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Section  1 
INTRODUCTION 

Satellite  imagery  suggests  that  large  portions  of  the  tropics  are  covered 
by  extensive  cirrus  cloud  systems  (Liou,  1986).  Tropical  cirrus  clouds  evolve 
during  the  life  cycle  of  the  mesoscale  convective  systems  and  are  modulated  by 
large-scale  disturbances.  Outflow  cirrus  clouds  from  tropical  cumulonimbi  appear 
to  be  maintained  in  a  convectively  active  state  by  radiative  flux  gradients 
within  the  clouds,  as  suggested  by  Danielson  (1982).  Extensive  anvils  are  likely 
to  become  radiatively  destabilized  by  cooling  at  tops  and  warming  at  bases .  This 
would  drive  convective  fluxes  which  in  turn  would  provide  an  upward  flvix  of  water 
vapor  within  the  cloud.  The  additional  moisture  at  cloud  top  levels  would 
promote  rapid  ice  crystal  growth  and  fallout.  Ackerman  et  al.  (1988)  have 
computed  radiative  heating  rates  in  typical  tropical  anvils.  The  heating  rate 
differences  between  the  cloud  bottom  and  top  ranges  from  30  to  200  K/day.  Lilly 
(1988)  has  analyzed  the  dynamic  mechanism  of  the  formation  of  cirrus  anvils  using 
a  mixed  layer  model,  and  has  shown  that  destabilization  of  the  layer  could  be 
produced  by  strong  radiative  heating  gradients.  The  importance  of  radiative 
processes  in  the  life  cycle  of  tropical  anvils  and  convective  systems  has  also 
been  illustrated  by  Chen  and  Cotton  (1988)  and  Dudhia  (1989). 

Clearly,  strong  radiative  heating  gradients  generated  by  tropical  cirrus 
anvils  would  have  a  significant  Impact  on  tropical  dynamic  and  thermodynamic 
processes  on  a  variety  of  scales  as  well  as  on  cloud  microphysical  processes. 
Those  processes  in  turn  modulate  the  cloud  evolution  and  the  associated  radiative 
heating  profiles.  Understanding  the  intricate  interactions  among  radiation, 
cloud  microphysics,  and  djmamics  requires  a  mesoscale  cloud  model  that  includes 
an  interactive  radiation  code.  We  are  in  the  process  of  merging  the  new  and 
comprehensive  radiation  parameterization  developed  by  Fu  and  Liou  (1992a, 


) 


Appendix),  specifically  designed  for  mesoscale  models,  with  the  cumulus  ensemble 
model  developed  by  Krueger  (1985  and  1988).  This  report  presents  some 
preliminary  results  from  the  integration  of  the  radiation  and  cloud  models. 

In  Section  2,  we  first  briefly  describe  the  cumulus  ensemble  model  (CEM) 
and  the  parameterization  of  three-phase  microphysics  processes.  Next,  we  outline 
the  parameterizations  of  scattering  and  absorption  processes  involving  absorbing 
gases  and  particulates  and  the  radiative  transfer  scheme .  In  Section  3 ,  we 
describe  the  initial  and  boundary  conditions  for  the  CEM  and  present  the 
microphysical  and  thermodynamic  fields  of  a  squall  line  system  simulated  by  the 
CEM.  Subsequently,  the  radiative  heating  rate  fields  and  cloud  radiative 
forcings  associated  with  this  squall  line  system  are  presented.  Finally,  a 
summary  is  given  in  Section  4. 
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Section  2 
MODEL  DESCRIPTION 

2.1  Cuffiulus  Ensemble  Model 

A  cunmlus  ensemble  model  (CEM)  Is  a  ntimerical  model  that  covers  a  large 
area  and  at  the  same  time  resolves  individual  cvimulus  clouds.  The  CEM  differs 
from  the  more  familiar  isolated  cloud  models  in  that  it  can  be  used  to  simulate 
the  response  of  a  cumulus  ensemble  to  a  prescribed  large-scale  condition.  It  can 
also  simulate  the  mesoscale  organization  of  convection.  According  to 
observations  (e.g.,  those  made  during  GATE),  a  three-dimensional  (3D)  CEM  is 
required  to  completely  describe  the  evolution  of  mesoscale  motions  as  well  as  the 
behavior  of  individual  cvimulus  clouds.  Since  a  tropical  anvil  cloud  associated 
with  a  mesoscale  convective  system  may  be  several  hundred  kilometers  wide  and 
last  for  several  hours,  the  CEM  should  cover  a  region  of  at  least  500  km  and 
should  be  Integrated  for  more  than  ten  hours  to  simulate  the  life  cycle  of  an 
anvil  complex.  Because  of  the  computational  effort  involved,  a  3D  CEM  simulation 
is  not  feasible  at  present,  especially  when  an  interactive  radiation  program  is 
included.  For  this  reason,  we  shall  confine  our  presentation  to  a  two- 
dimensional  (2D)  CEM.  Many  pioneering  studies  have  demonstrated  that  a 
significant  physical  understanding  regarding  the  mesoscale  convective  system  can 
be  derived  from  simulations  using  a  2D  CEM  model  (Llpps  and  Hemler,  1986;  Tao  et 
al.  1987;  Rotunno  et  al.  1988;  Krueger,  1988;  Xu  and  Krueger,  1991). 

In  the  present  study  we  use  the  2D  CEM  developed  by  Krueger  (1988)  .  In  the 
following,  we  present  the  dynamic  structure  of  the  model,  which  is  based  on  the 
anelastic  set  of  equations.  In  Cartesian  coordinates  for  the  case  of  a  slab 
symmetry,  the  basic  equations  for  momentvim,  continuity,  thermodynaimics ,  and  water 
mixing  ratios  may  be  written  as  follows: 
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^<“>  -  »  f<V>, 


(2.1) 


^  <V>  -  -^<u"v">-^  T^(Pn<V"w">)-f  [<U>-Il(z)  ]  , 
^  Po  dz  “ 


■at 


(2.2) 


<«1> 


+  0.61<qyl>-<q^>-<qp— <qi>-<qs>-<qg> 


(2.3) 


^(,0<U»  .  .^(P.<»»  ■  0, 


^<S>  -  -.l^(Po<w''fl''»-.^<u-S''>-w(z,t).^<tf>  ♦  s,  +  R,. 


Po 

-w(z,t).^<q^  ,  i>  +  ,,, 


i  > 


•3r£(^0'''^r.s,g<‘lr,8,g>)  Sr,»,g> 


where  the  potential  temperature  is  given  by 


0  .  T 


the  virtual  potential  temperature  is 


(C' 


the  dimensionless  pressure  is 


6^  -  0(1  ♦  0.61q^), 


’  ■  i^r 


(2.4) 


(2.5) 


(2.6) 


(2.7) 


and  the  hydrostatic  equation  is 
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^”’o  ,  _  g 

CpffyQ 

In  the  preceding  equations,  the  angle  brackets  denote  the  ensemble  mean 
(equivalent  to  an  average  over  the  y  -  coordinate ) ;  the  double  primes  represent  the 
departure  from  the  ensemble  mean;  the  variables  with  subscript  0  refer  to  the 
reference  state,  which  is  hydrostatic  and  is  a  function  of  the  height,  z,  only; 
the  subscript  1  denotes  the  departure  from  the  reference  state;  u,  v,  and  w  are 
the  three-dimensional  velocity  components;  t  is  the  time;  x  and  z  are  the 
horizontal  and  vertical  coordinates;  Cp  is  the  specific  heat  at  constant 
pressure;  f  is  the  Coriolis  parameter;  g  is  the  acceleration  of  gravity;  R  is  the 
gas  constant  for  dry  air;  p  is  the  air  density;  T  is  the  air  temperature;  p  is 
the  pressure;  poo  is  a  constant  reference  pressure;  w  is  the  prescribed  large- 
scale  vertical  velocity;  u,  is  the  prescribed  geostrophic  wind;  qv,c,i,r.s,g  *^he 
mixing  ratios  and  ^  ^  j.  ,  g  represent  the  corresponding  sources/sinks  due  to 
water  phase  changes,  in  which  the  subscripts  v,  c,  i,  r,  s,  and  g  denote  water 
vapor,  cloud  water,  cloud  ice,  rain,  snow,  and  graupel,  respectively;  Sg  and  Rj 
represent  the  sources/sinks  of  heat  due  to  water  phase  transitions  and  radiative 
processes,  respectively,  as  expressed  by  changes  in  potential  temperature;  and 
Vj.  ,  g  are  the  mass-weighted  mean  terminal  velocities  of  rain,  snow  or  graupel 
(Strivastava,  1967).  In  the  preceding  equations,  the  effect  of  the  large-scale 
vertical  velocity  is  included  only  in  the  equations  for  potential  temperature  and 
suspended  water.  Although  the  direct  effects  of  the  large-scale  vertical 
velocity  on  cumulus  developments  are  negligible,  it  will  influence  cumulus 
convection  through  the  temperature  and  water  vapor  fields. 

Nineteen  turbulent  flux  terms  are  required  to  close  Eqs .  (2.1) -(2. 7).  The 
turbulent  fluxes  for  rain,  snow,  and  graupel  are  set  to  zero,  while  <u"v">  and 
<w"v">  are  determined  by  using  the  Mellor-Yamada  (1974)  level  one  equation  as 
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presented  by  Xu  and  Krueger  (1991) .  The  remaining  turbulent  fluxes  are  predicted 
using  third  moment  closure  (Krueger,  1988).  In  the  turbulence  equations,  the 
turbulent  fluxes  of  the  mixing  ratio  of  suspended  water,  qw(“qv  +  ^Ic  ^  9i). 
calculated;  fluxes  of  the  individual  components  are  not  required.  Similarly, 
turbulent  fluxes  of  the  liquid-ice  water  static  energy,  S,  -  CpT  +  L^qy  -  Lfq^  - 
Lcqc  +  gz,  where  and  Lf  are  the  condensation  and  freezing  latent  heats,  are 
calculated;  fluxes  of  6  are  not  required.  The  third-moment  closure  has  a 
distinct  advantage  over  simpler  closures  because  it  is  more  general  and  should 
therefore  have  a  wider  range  of  validity.  In  particular,  it  can  be  utilized  to 
improve  the  simulation  of  boundary -layer  turbulence  and  the  treatment  of  in-cloud 
turbulence . 

2.2  Mlcrophyslcal  Parameterization 

The  present  model  utilizes  a  bulk  parameterization  to  represent  cloud 
microphysical  processes,  in  which  the  hydrometeors  are  categorized  into  five 
types:  cloud  droplets ,  ice  crystals,  raindrops,  snow  (or  aggregates) ,  and  graupel 
(or  hail) .  The  cloud  microphysics  parameterizations  follow  the  procedures 
developed  by  Lin  et  al.  (1983)  and  Lord  et  al.  (1984)  and  are  briefly  described 
in  the  following.  The  microphysics  code  was  provided  by  Dr.  Stephen  Lord. 

It  is  convenient  to  write  the  source-sink  term  in  Eqs.  (2.6)  and  (2.7)  as 

S,  -C,  ^P,,  (2.8) 

where  the  subscript  x  denotes  v,  c,  i,  r,  s,  or  g;  C  is  the  source-sink  term  due 
to  the  various  possible  interactions  among  water  vapor,  cloud  water,  and  cloud 
ice  (Cr  "  Cj  —  Cg  —  0) ;  and  P  is  the  production  terms  due  to  a  large  number  of 
individual  microphysical  processes,  which  are  related  to  precipitating  particles 
and  are  represented  by  using  bulk  parameterizations.  These  microphysical 


Table  1. 


Microphysical  processes  represented  by  bulk  parame ter izat ions  (after 
Lord  et  al . ,  1984) . 


Symbol 

Meaning 

R*vp 

Evaporation  of  rain 

^-aut 

Autoconversion  (collision-coalescence)  of  cloud  water  to 

form  rain 

^acw 

Accretion  of  cloud  water  by  rain 

^acB 

Accretion  of  snow  by  rain 

Ssub 

Sublimation  of  Snow 

Sdap 

Depositional  growth  of  snow 

^aut 

Autoconversion  (aggregation)  of  cloud  ice  to  form  snow 

^Bllt 

Melting  of  snow 

Sfi 

Transfer  rate  of  cloud  ice  to  snow  through  depositional 
growth  of  Bergeron  process  embryos 

Sf 

Transfer  rate  of  cloud  water  to  snow  through  Bergeron 
process  (deposition  and  riming) 

^acw 

Accretion  of  cloud  water  by  snow 

Saei 

Accretion  of  cloud  ice  by  snow 

^aer 

Accretion  of  rain  by  snow 

Gfub 

Sublimation  of  graupel 

®aut 

Autoconversion  (collision-aggregation)  of  snow  to  form 
graupel 

Gfr 

Probabilistic  freezing  of  rain  to  form  graupel 

Ginlt 

Melting  of  graupel  to  form  rain 

Gnat 

Wet  growth  of  graupel 

Gacw 

Accretion  of  cloud  water  by  graupel 

Gael 

Accretion  of  cloud  ice  by  graupel 

G,cr 

Accretion  of  rain  by  graupel 

Gaea 

Accretion  of  snow  by  graupel 

processes  are  described  in  Table  1.  The  bulk  representation  of  the  individual 
production  terms  can  be  found  in  Lin  et  al.  (1983).  The  total  production  term 
for  six  types  of  water  substance  may  be  written  as  follows; 

(i)  If  the  temperature  is  below  0®C  we  have 
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*  -(l-«i)R.vp-(l-«i)S.ub-5iSd,p-(l-fi:)G.^,  (2.9a) 

Rc  “  “R«ut“R«cw"^fw“^acw“^acw>  (2.9b) 

Pi  *  ,  (2.9c) 

P,  =  (l-Si)R.^  +  R.„t  ♦  R.cw-S.cr-G£r-G.cr.  (2.9d) 

Pa  “  “(^~^2)Rae8  ■*■  ^2^acr  "*■  ■*■  ^l^dap  (2.9e) 

*  ^aut  *  ^£w  *  ^tl  *  ^acw  '*'  ^acl  “Gaut~Gaca  > 


P,  =  (l-i2)R.e.  ♦  (l-«2)Sacr  ^  +  (l-5l)G.^  * 

+  inin(G^jj^ ,  Gjj.^  +  G,(,j^  •♦•  G^g^  +  G^^,)  . 

(ii)  If  the  temperature  is  above  0*C  we  have 


P^  - -(l-«i)R.^.  (2.10a) 

Pq  *  ~Raut”R'acw'"^aew~G,cw>  (2.10b) 

Pi  -  0.  (2.10c) 

Rr  '  (l~^l)Revp  ■*■  Raut  ■*■  Raew  *  ^aca  ~  ^alt  *  ^acw  ~  G^it  "*■  ^acw>  (2.10d) 

R.  •  S„i,-G,„-R.„,  (2.10e) 

Rs  =G„i,  ^G.„.  (2.10f) 

In  Eqs.  (2.9a)  -  (2.9f),  5^  -  1  within  cloud  (q^  +  q^  >  0)  ,  5^  -  0  otherwise;  S2 
-  1  for  q^  and  q,  <  0.1  g/kg,  and  S2  ~  0  otherwise.  In  Eqs.  (2.10a)  -  (2.10f) 
S I  *  S 2  *  0. 

The  determination  of  C  terms  follows  Lord  et  al.  (1984).  The  following 


basic  assumptions  are  made.  First,  all  supersaturated  vapor  condenses/deposits, 
and  is  partitioned  betT<‘.en  cloud  water  and  cloud  ice  as  a  linear  function  of 
temperature.  Second,  no  ice  is  produced  at  T  >  O'C,  while  no  cloud  water  is 
produced  for  T  <  -40“C,  at  which  homogeneous  nucleation  is  assumed  to  occur. 
Third,  the  saturation  water  vapor  mixing  ratio  is  delined  as 
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for  Qc  +  qi  >  0 
for  qc  =  Qi  =  0 , 


(2.11) 


where  q,*  and  q^*  are  saturation  water  vapor  mixing  ratios  with  respect  to  liquid 
water  and  ice,  respectively.  Fourth,  the  supersaturation  is  removed  at  the  end 
of  each  time  step;  this  is  referred  to  as  the  saturation  adjustment  scheme. 
Fifth,  the  ice  crystal  nucleation  process,  Ijh,,  as  described  by  Lin  et  al.  (1983) 
is  retained  as  one  of  the  individual  microphysical  processes  to  account  for  the 
conversion  of  cloud  water  to  cloud  ice  through  the  Bergeron  process  before  the 
saturation  adjustment. 

Mass-weighted  mean  terminal  velocities  for  rain  (Vj.) ,  snow  (Vg) ,  and 
graupel  (Vg)  are  derived  by  Lin  et  al.  (1983).  Typical  values  of  these 
velocities  are  5,  2,  and  10  ms"^  for  rain,  snow,  and  graupel,  respectively.  The 
$9  term  in  the  thermodynamic  equation  can  be  written  as 


Sj  =  Pj+C(,,  (2.12) 
where  Pg  is  the  potential  temperature  change  due  to  the  microphysical  processes 
listed  in  Table  1  and  Cj  is  due  to  the  phase  changes  among  water  vapor,  cloud 
droplets ,  and  cloud  ice .  We  have 


^  [(l-5i)4R,^  >  L£(S„i,  +  -  R.c.)]  for  T  >  0-C, 

"  ^(Sr«  *  S.„  G,,  (2.13) 

[  (^ “  ^sub  "*■  ^l^dap 

+  (l-5i)Gg^]}  forT<0*C, 


where  k  -  R/Cp,  and  L^,  L^,  and  L,  are  the  latent  heats  of  condensation,  freezing, 
and  sublimation,  respectively.  Cg  can  be  expressed  as 
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for  T  >  CC 


1  fPoor  I^^qc  "  Mqi  forT>0»C 

^  I  p  J 

^  £<.rT<0-C, 


(2.14) 


where  Aq^  and  Aq^  are  the  production  of  cloud  water  and  ice  calculated  by  the 
saturation  adjustment  scheme  and  At  the  time  step. 


2.3  Radiative  Transfer  Scheme 


The  Rg  term  in  Eq.  (2.5)  can  be  written  as 


(2.15) 


where  (3T/dt)R  represents  the  radiative  heating  rate,  which  can  be  computed  from 


fail  .  _  1  d(F^-F^) 

[3t^  ”7^  cE  ’ 


(2.16) 


with  F’  and  F*  the  net  upward  and  downward  fluxes  covering  both  solar  and  thermal 
infrared  spectra.  Radiative  fields  in  the  cloud  are  strongly  modulated  by  the 
cloud  microphysical  structure.  We  have  developed  a  numerically  stable  and 
computationally  efficient  scheme  for  the  calculation  of  upward  and  downward 
fluxes.  This  radiation  scheme  is  formulated  in  such  a  manner  that  absorption  and 
scattering  by  both  molecules  and  hydrometeors  can  be  treated  consistently. 


2.3.1  Parameterization  of  the  Single-Scattering  Properties  of  Hydrometeors 

The  single-scattering  properties,  including  the  extinction  coefficient,  )9, 
the  single  scattering  albedo,  w,  and  the  scattering  phase  function,  P(cos0)  ,  must 
be  known  before  radiative  transfer  calculations  can  be  performed.  The  single¬ 
scattering  properties  of  hydrometeors  are  related  to  the  thermodynamic  phase  (ice 
or  water),  shape,  and  size  distribution. 
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a.  Ice  Crystals 

In  this  study,  ice  crystals  are  assumed  to  be  hexagonal  columns  and  plates 
with  length,  L,  and  width,  D.  To  the  extent  that  the  scattering  of  light  is 
proportional  to  the  cross  section  area  of  randomly  oriented  hexagonal  ice 
crystals,  we  may  define  a  mean  effective  size  to  represent  ice  crystal  size 
distribution,  n(L) ,  in  the  form 

D,  *  J  D-DL  n(L)  dL/J  DL  n(L)  dL.  (2.17) 

Based  on  physical  principles,  as  discussed  by  Fu  and  Liou  (1992a),  the 
extinction  coefficient  and  single -scattering  albedo  for  ice  crystals  can  be 
parameterized  by  the  following: 

^-IWcfa„/D;.  (2.18) 

nHo 

3 

“  *  E  »  (2.19) 

h»0 

where  IWC  is  the  ice  water  content,  and  a^  and  b„  are  certain  coefficients.  The 

phase  function  is  usually  expanded  in  temiis  of  Legendre  polynomials  P,  in  the 
M  _ 

form  P(cos0)  -  P|  (cos0) ,  where  the  expansion  coefficients,  w, ,  can  be 
expressed  by 


(l-f,)w^  +  fj(2i  +  l)  for  solar  ^2  20) 

W,  »  *  \  •  J 

(2i+l)g'  for  IR, 

where  «,*  represent  the  expansion  coefficients  without  the  incorporation  of  the 
5-function  transmission  through  parallel  planes  at  0  -  0  (Takano  and  Liou,  1989), 
f,  denotes  the  forward  contribution  due  to  the  5 -function  transmission,  and  g  is 
the  asymmetry  factor,  w,*,  f^,  and  g  can  be  parameterized  by 
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Table  2.  Characteristics  of  the  11  ice  crystal  size  distributions. 


Particle  Size 
Distribution 

Ice  Water 

Content  (g  m'^) 

Mean  Effective  Size 
(pm) 

Cs* 

4.765  e-3 

41.5 

Ci*  Uncinus 

1.116  e-1 

123.6 

Ci**  (cold) 

1.110  e-3 

23.9 

Ci**  (warm) 

9.240  e-3 

47.6 

Ci**  (T— 20“C) 

8.613  e-3 

57.9 

Ci**  (T— 40-0 

9.177  e-3 

64.1 

Ci**  (T— 60“C) 

6.598  e-4 

30.4 

Ci***  (Oct.  22) 

1.609  e-2 

104.1 

Ci***  (Oct.  25) 

2.923  e-2 

110.4 

Ci***  (Nov.  1) 

4.968  e-3 

75.1 

Ci***  (Nov.  2) 

1.406  e-2 

93.0 

♦Heymsfield  (1975);  **Heymsfield  and  Platt  (1984);  ***FIRE  (1986). 


Ut 


3 

E 

n«0 


n«0 


(2.21)  • 


g(IR) 


3 


D." 


n*0 


(2.22) 


where  c^  , ,  d„,  and  c'„  are  certain  coefficients.  In  the  thermal  infrared 
wavelengths,  it  suffices  to  use  the  asymmetry  factor  via  the  Henyey- Greens te in 
function  to  represent  the  phase  function  because  the  halo  and  6 -transmission  peak 
features  in  the  P(cos0)  are  largely  suppressed  due  to  strong  absorption. 

The  extinction  coefficient,  single-scattering  albedo,  and  expansion 
coefficients  of  the  phase  function  for  the  11  observed  ice  crystal  size 
distributions  (Table  2)  have  been  computed  from  a  geometric  ray-tracing  method 
for  hexagonal  ice  crystals  (size  parameter  >  30)  and  from  a  Mie-type  solution  for 
spheroids  (size  parameter  <  30).  The  coefficients,  a^,  b^,  c'„,  c„  ,  and  dj,  are 
obtained  by  numerical  fitting  to  the  data  computed  from  the  "exact"  computations . 
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Table  3.  Characteristics  of  the  eight  droplet  size  distributions. 


Cloud  Type* 

Liquid  Water 

Content  (g  m'^) 

Mean  Effective  Radius 
ipm) 

St  I 

0.22 

5.89 

St  II 

0.05 

4.18 

Sc  I 

0.14 

5.36 

Sc  II 

0.47 

9,84 

Ns 

0.50 

9.27 

As 

0.28 

6.16 

Cu 

1.00 

12.10 

Cb 

2.50 

31.23 

^Stephens  (1978) 


The  preceding  parameterlzatlons  have  relative  accuracies  within  -1%. 

For  radiative  transfer  calculations,  nonspherical  ice  crystals  have  been 
frequently  approximated  by  ice  spheres  with  equivalent  areas  (see  e.g.  Stackhouse 
and  Stephens,  1991).  Equivalent- area  spheres  scatter  more  light  in  forward 
directions  and  have  smaller  single-scattering  albedos  than  nonspherical  ice 
crystals.  As  a  result,  the  assumption  that  ice  crystals  are  spheres  leads  to  a 
significant  underestimation  for  the  solar  albedo  of  cirrus  clouds.  As  shown  by 
Fu  and  Llou  (1992a),  the  present  parameterization  can  be  used  to  reasonably 
interpret  the  observed  IR  emisslvities  and  solar  albedo  involving  ice  clouds. 
b.  Water  Droplets 

For  water  clouds,  we  may  define  a  mean  effective  radius,  r, ,  to  represent 
the  water  droplet  size  distribution,  n(r),  with  respect  to  radiative  calculations 
as  follows: 


r,  »  r  r^  n(r)dr/r  r^  n(r)dr. 


(2.23) 
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The  single-scattering  properties  of  water  droplets  can  be  calculated  exactly  by 
using  the  Mie  theory.  Eight  water  cloud  types  presented  by  Stephens  (1978)  are 
used  in  the  present  study  and  are  listed  in  Table  3.  By  using  the  single¬ 
scattering  properties  of  the  eight  water  clouds,  the  following  parameterizations 
have  been  developed  (Fu,  1991): 


LWC 


/3i  ^  ^2/LWC2  -  ^i/LWCi  ^ 

'  Y 

1  _  1 

-  1/r-  -  ^ 

J. 

(2.24) 


W  *  Wi  ♦ 


^.2  - 


(2.25) 


Si*  -fi—  (r.  -  r,^). 

•2 


(2.26) 


where  LWC  is  the  liquid  water  content.  The  terms  wj,  gi)  and  (y32,  W2>  S2) 

denote  the  single -scattering  properties  for  two  cloud  types  listed  in  Table  3 
with  (LWCi,  r,i)  and  (LWC2,  r,2) .  respectively,  so  that  r,i  and  r,2  are  the  pair 
of  mean  effective  radii  closest  to  r,  and  r,i  <  r,  <  r,2.  The  scattering  phase 
functions  for  water  clouds  are  approximated  by  the  Henyey- Greens te in  function 
through  the  asymmetry  factor, 
c.  Raindrops,  Snow,  and  Craupel 

The  size  distributions  for  these  hydrometeors  used  in  the  bulk 
parameterization  are  in  the  forms 


n,(D)  -  no,  exp(-A,D,)  ,  (2.27) 

where  the  subscript  x  denotes  r,  s,  or  g;  no  is  the  intercept  parameter;  A  is  the 
slope  parameter;  and  D  is  the  diameter  of  the  hydrometeors.  The  intercept 
parameters,  as  given  by  Lord  et  al.  (1984),  are  no,  -  0.22,  nQ,  -  0.03,  and  nog 
-  4xl0'*cm"*.  The  slope  parameter  is  given  by 
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*  pq* 


where  -  1,  p,  -  0.1,  and  pg  -  0.3  g/cm^;  and  pq^  is  the  rain,  snow,  or  graupel 
water  content  (g/cm®)  .  Setting  pq^  -  0.5x10'®  g/cm®,  we  have  -  34.3,  A,  -  11.7, 
and  Ag  -  5.24  cm'^.  The  corresponding  mean  effective  radii  from  Eq.  (2.27)  are 
0.044,  0.128,  and  0.286  cm  for  rain,  snow,  and  graupel,  respectively.  These  mean 
effective  radii  are  determined  from 


The  single -scattering  properties  u>\  and  g',  for  rain,  snow,  and  graupel 
with  the  size  distributions  defined  in  Eq.  (2.27)  for  the  given  A^  above  are 
calculated  by  the  Mie  theory.  For  a  water  content,  pqx,  from  the  cloud  model, 
we  have 


Pqr^r 

1.72x10'= 


(2.28) 


pq.^s 

0.585x10 


> 


(2.29) 


a- 

0.262x10'=  *  ’ 


(2.30) 


where  pq,  and  A,  are  in  units  of  g  cm'=  and  cm'^,  respectively.  The  single- 
scattering  albedo  and  asymmetry  factor  are  assumed  to  be  equal  to  w’  and  g' 
computed  from  Mie  calculations  regardless  of  the  size  distribution  dependence. 


2.3.2  Parameterization  of  Nongray  Gaseous  Absorption 

The  solar  and  thermal  infrared  spectra  are  divided  into  a  number  of  bands 
depicted  in  Table  4.  Mean  values  for  the  scattering  properties  of  hydrometeors, 
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Table  4.  Spectral  division  used  in  the  parameterization. 


Solar  Spectrum 

Infrared  Spectrum 

Band 

(i)  Central  A 

Band  Limit 

Band 

(i) 

Central  A 

Band  Limit 

(pm) 

(pm) 

(pm) 

(cm- 1) 

1 

0.55 

0.2-0. 7 

7 

4.9 

2200-1900 

2 

1.0 

0.7-1. 3 

8 

5.6 

1900-1700 

3 

1.6 

1. 3-1.9 

9 

6.5 

1700-1400 

4 

2.2 

1.9-2. 5 

10 

7.6 

1400-1250 

5 

3.0 

2. 5-3. 5 

11 

8.5 

1250-1100 

6 

3.7 

3. 5-4.0 

12 

9.6 

1100-980 

13 

11.3 

980-800 

14 

13.7 

800-670 

15 

16.6 

670-540 

16 

21.5 

540-400 

17 

30.0 

400-280 

18 

70.0 

280-10 

solar 

irradiances,  and 

Planck  functions 

are  used  in 

each  band. 

In  the  solar 

spectrum,  absorption  due  to  H2O  (2500-14500  cm"^) ,  O3  (in  the  ultraviolet  and 
visible),  CO2  (2850-5250  cm*^)  ,  and  O2  (A,  B,  and  7  bands)  is  accounted  for  in  the 
radiation  scheme.  In  the  infrared  spectrum,  we  include  the  major  absorption 
bands  of  H2O  (0-2200  cm'^)  ,  CO2  (540-800  cm’^),  O3  (980-1100  cm’^)  ,  CH^  (1100-1400 
cm'^)  ,  and  N2O  (1100-1400  cm'^)  .  The  continuum  absorption  of  H2O  is  incorporated 
in  the  spectral  region  280-1250  cm'^. 

Nongray  gaseous  absorption  is  parameterized  based  on  the  correlated  k- 
distribution  method  developed  by  Fu  and  Liou  (1992b).  In  this  method,  the 
cumulative  probability,  g,  of  the  absorption  coefficient,  k„,  in  a  spectral 
interval,  Ai/,  is  used  to  replace  the  frequency,  1/,  as  an  independent  variable. 
This  leads  to  an  immense  numerical  simplification,  in  which  about  ten  thousand 
frequency  intervals  can  be  replaced  by  a  few  g  intervals.  Using  a  minimum  number 
of  g  intervals  to  represent  the  gaseous  absorption  and  to  treat  overlap  within 


each  spectral  interval,  121  spectral  calculations  are  required  for  each  vertical 
profile.  Compared  with  results  from  a  LBL  program,  the  parameterizations  achieve 
an  accuracy  within  0.1  k/day  for  heating  rates  and  0.5%  for  fluxes. 

2.3.3  Parameterization  of  Radiative  Transfer 

For  parameterizations  of  radiative  fluxes,  we  use  the  5 -four- stream 
approximation  developed  by  Liou  et  al.  (1988).  For  a  homogeneous  layer,  an 
analytic  solution  can  be  derived  explicitly  for  this  approximation  so  that  the 
computer  time  Involved  is  minimal.  In  order  to  apply  this  approach  to  the 
thermal  infrared  radiative  transfer,  the  Planck  function  is  expressed  in  terms 
of  optical  depth,  r,  in  the  form  a»exp(br),  where  a  and  b  are  coefficients 
determined  from  the  top-  and  bottom-layer  temperatures.  Since  the  direct  solar 
radiation  source  also  has  exponential  function  form  in  terms  of  optical  depth, 
the  solution  of  the  j- four-stream  approximation  for  IR  wavelengths  is  the  same 
as  that  for  solar  wavelengths.  For  application  to  a  nonhomogeneous  atmosphere 
we  divide  this  atmosphere  into  N  layers  within  which  the  5 -four- stream  scheme  can 
be  applied.  The  unknown  coefficients  in  the  analytic  solution  for  the  radiative 
transfer  equation  are  determined  following  the  procedure  described  in  Liou 
(1975).  The  total  single -scattering  properties  due  to  the  combined  contributions 
of  Rayleigh  scattering,  nongray  gaseous  absorption,  and  scattering  and  absorption 
by  hydrometeors  can  be  evaluated  following  the  procedures  described  by  Fu  and 
Liou  (1992a). 

2.3.4  Comparison  with  ICRCCM  Results 

Results  computed  from  the  present  radiation  scheme  are  compared  with  those 
presented  from  the  Intercomparison  of  Radiation  Codes  in  Climate  Models  (ICRCCM) 
program  for  cloudy  conditions.  In  the  ICRCCM,  six  sets  of  radiation  calculations 
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Table  5.  Statistics*  on  the  downward  IR  surface  flux  calculations  under 

overcast  cloudy  conditions.  Results  computed  from  the  present  # 

radiation  scheme  are  depicted  in  the  parentheses . 


Cloud  Type 

Top 

Height 

(km) 

LWP 

(g 

Number 
of  Model 

Median 

^  surface 

(W  m'^) 

Range 

(%) 

rms  Diff. 

(%) 

41 

CS 

2 

10 

10 

399 

(395) 

7.0 

2.1 

CL 

2 

10 

7 

387 

(368) 

9.0 

2.9 

CL 

2 

200 

14 

413 

(411) 

1.9 

0.6 

9 

CS 

13 

10 

11 

360 

(361) 

6.4 

1.9 

CL 

13 

10 

6 

358 

(353) 

5.6 

2.0 

• 

CL 

13 

200 

13 

361 

(362) 

6.7 

2.0 

♦After  Ellingson  and  Fouquart  (1990).  ^ 


were  performed  for  overcast  conditions  with  the  aim  of  testing  the  sensitivity 

of  the  radiation  program  to  the  drop  size  distribution,  the  location  of  the  cloud  # 

top,  and  the  cloud  LWC.  The  cloud  thickness  was  assumed  to  be  1  km,  with  the 

cloud  tops  set  at  2  and  13  km  for  low  and  high  clouds,  respectively.  Two  droplet 

size  distributions  were  selected;  one  with  small  droplets  (CS,  r,  -  5.36  /xm) ,  # 

while  the  other  with  large  droplets  (CL,  r,  -  31.23  urn).  The  LWC  was  specified 

to  be  either  10  g  m"^  (nonblack  clouds  for  both  CS  and  CL)  or  200  gm'^  (near-black 

clouds  for  CL)  .  The  calculations  were  performed  with  one  cloud  layer  present  in  9 

a  midlatitude  atmosphere.  For  solar  radiation  computations,  the  solar  zenith 

angle  and  the  surface  albedo  are  set  at  30®  and  0.2,  respectively. 


Table  5  presents  the  statistics  determined  from  the  ICRCCM  for  downward  IR 
surface  fluxes.  Results  from  the  present  radiation  scheme  are  included  in  the 
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parentheses.  For  atmospheres  containing  near-black  clouds,  the  fluxes  computed 
from  the  present  scheme  agree  with  the  median  values  provided  by  the  ICRCCM 
within  2  W/m^.  The  present  results  also  agree  well  with  those  of  ICRCCM  in  the 
case  of  CS.  Differences  of  up  to  19  W/m"^  in  the  surface  flux  are  seen  for  the 
case  of  CL  with  IWP  of  10  g  m"^  and  cloud  top  of  2  km.  The  downward  surface  flux 
from  the  ICRCCM  is  not  sensitive  to  r,  because  some  models  have  not  explicitly 
considered  the  dependence  of  the  cloud  radiative  properties  on  r^  (see,  e.g.  Liou 
and  Wittman,  1979;  Stephens,  1978). 

Table  6  shows  the  results  for  solar  radiation.  The  rms  differences  and 
total  ranges  from  ICRCCM  are  quite  large,  indicating  the  difficulty  of  adequately 
modeling  the  effects  of  multiple  scattering.  Differences  between  the  present 
results  and  the  medians  from  ICRCCM  are  much  smaller  than  the  rms  differences , 
especially  for  total  atmospheric  absorption. 


# 
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Section  3 

SIMULATION  OF  TROPICAL  CONVECTION 
3.1  Botmdary  and  Initial  Conditions 

In  the  CEM  simulation  analyzed  here,  we  use  a  horizontal  domain  of  1024  km 
and  a  horizontal  grid  size  of  2  km.  The  vertical  domain  used  is  19  km  with  a 
stretched  grid  consisting  of  33  layers.  Near  the  surface,  the  grid  interval  is 
100  m,  while  it  is  1  km  near  the  model  top.  The  lateral  boundary  condition  is 
cyclic  and  the  upper  and  lower  boundaries  are  rigid.  Numerical  simulations  are 
carried  out  over  the  ocean.  The  sea  surface  temperature  is  fixed  at  299.9  K. 
All  surface  turbulent  fluxes  are  diagnosed  by  using  the  flux-profile 
relationships  given  by  Deardorff  (1972). 

In  Eqs.  (2.5)  and  (2.6),  the  terms  that  include  the  large-scale  vertical 
velocity  are  prescribed.  These  terms  are  horizontally  uniform.  They  vary  with 
height  and  time  according  to 

f(z.t)  -  f(z)(l  +  cos(2xt/T)]/2, 

where  T  (-  27  h)  is  the  period,  and  f(z) ,  representing  typical  GATE  phase-IIl 
mean  profiles,  is  described  in  Xu  et  al.  (1992).  The  time- independent  x- 
component  of  the  geostrophic  wind  is  prescribed;  it  is  identical  to  that  used  by 
Xu  and  Krueger  (1991)  with  shear.  The  sheared  profile  of  the  geostrophic  wind 
is  typical  of  the  11  September  1974  squall  line  environment  observed  during  the 
GATE  phase -II I. 

The  initial  thermodynamic  conditions  are  horizontally  uniform.  Cloud 
fields  are  initiated  by  introducing  small,  random  temperature  perturbations  into 
the  lowest  model  layer  after  the  first  30  minutes  of  integration.  The  initial 
thermodynamic  state  used  in  this  simulation  is  identical  to  that  used  in  Xu  and 


Krueger  (1991).  The  numerical  simulation,  with  a  Coriolis  parameter  of  15“N,  was 
run  for  11  days  with  a  time  step  of  10  seconds. 

3.2  Thermodynamic  and  Cloud  Microphysical  Fields 

Figure  1  shows  the  time  evolution  of  the  cloud  top  temperature  from  3  to 
7  days,  indicated  by  a  linear  gray  scale;  white  represents  200  k,  while  black 
denotes  300  k.  Cirrus  anvils  associated  with  cumulonimbi  appear  white.  From 
Fig.  1  we  see  that  cumulus  convection  is  organized  into  long-lived  mesoscale 
systems.  These  systems  have  convective  bands  of  small  horizontal  extent  that  are 
accompanied  by  cirrus  anvils  behind  the  bands.  Some  midlevel  stratiform  clouds 
also  appear  behind  the  bands.  To  study  the  effects  of  the  tropical  mesoscale 
convective  system  on  the  radiative  fields,  the  simulated  thermodynamic  and  bulk 
hydrometeor  fields  are  obtained  from  a  6-hour  average  from  4.75  to  5  days.  In 
the  next  section  we  will  present  the  radiation  budgets  of  the  mesoscale 
convective  system  computed  from  these  fields.  The  6-hour  average  thermodynamic 
and  microphysical  fields  are  presented  below. 

Figure  2  shows  the  detailed  cross  sections  (x-z)  of  the  simulated  cloud  and 
precipitation  fields.  The  leading  convective -stratiform  area  has  a  large 
horizontal  extent  of  -320  km.  The  anvil  cloud  associated  with  the  mesoscale 
convective  system  has  a  horizontal  scale  of  -250  km,  which  extends  from  the 
middle-front  to  the  back  of  the  system.  The  system  moves  westward  (from  right 
to  left)  and  consists  of  significant  precipitation  that  covers  a  region  of  -100 
km.  To  examine  the  spatial  distribution  of  the  radiation  field  associated  with 
a  mesoscale  convective  system,  we  shall  confine  our  presentation  to  the 
horizontal  domain  from  120  to  620  km. 

Figures  3,  4,  and  5  show  the  x-z  sections  of  the  ice  phase  mixing  ratio, 
water  phase  mixing  ratio,  and  total  hydrometeor  mixing  ratio,  respectively.  The 
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Fiq.  1.  Hovmbllec  diagram  (x-t  section)  of  the  cloud  top  temperature,  based  on 
15-minute  average.  Cloud  top  temperature  is  denoted  by  a  linear  gray 
scale;  white  and  black  represent  200  and  300  K,  respectively. 
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section  (x-z)  of  the  mixing  ratio  for  total  ice  species  (cloud  ice,  snow,  and  graupel).  The  contours 
logarithmic  scale  with  an  interval  of  10®.  The  lowest  contour  corresponds  to  10'^. 
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anvil  base  (Fig.  3)  is  estimated  to  be  at  about  3-7  km.  The  anvil  base  is  not 
fixed  because  it  consists  of  melting  snow/graupel .  The  anvil  as  shown  in  Fig. 
3  is  extremely  thick,  -4-12  km,  with  cloud  tops  located  at  -10-15  km.  The 
maximum  total  ice  water  mixing  ratio  is  -1  g/kg  in  this  case.  The  maximum  liquid 
water  mixing  ratio  (Fig.  4)  is  -2.4  g/kg,  which  is  largely  due  to  raindrops. 
Heavy  rain  occurs  just  behind  the  gust  front  (the  left  side  of  the  deep 
cumulonimbi  which  form  the  squall  line  is  the  location  where  the  principal 
updraft  occurs),  as  seen  from  Fig.  5.  The  deep  cumulonimbi  continually  propagate 
into  the  ambient  air  ahead  where  new  growth  occurs,  while  older  towers 
successively  Join  the  main  anvil  mass  (Zipser,  1977).  The  rain  observed  in  the 
310-350  km  zone  comes  from  the  anvil  base,  where  very  few  low  clouds  exist  (See 
Fig,  2  and  5),  Figure  6  shows  the  vertical  distribution  of  the  mixing  ratios  for 
total  liquid  water,  ice  water,  and  hydrometeor  at  217,  317,  and  417  km.  The 
microphysical  distributions  at  those  three  positions  represent  the  pictures  in 
front  of  the  squall  line,  at  the  center  of  the  cloud  area,  and  in  the  rear  of  the 
system,  respectively.  At  217  km,  the  total  hydrometeor  mixing  ratio  is  dominated 
by  water  phase,  while  at  417  km,  ice  phase  prevails.  At  x  -  317  km,  the  total 
liquid  water  is  largely  due  to  raindrops  which  fall  from  the  anvil.  The  total 
liquid  water  mixing  ratios  at  217  and  417  km  are  largely  due  to  cloud  drops. 
From  Fig.  6,  it  is  seen  that  the  mixing  ratio  of  ice  phase  reaches  its  maximum 
in  the  middle  of  the  anvil,  while  water  phase  shows  a  maximum  mixing  ratio  near 
the  cloud  top. 

Figure  7  shows  the  first  model  layer  temperature  above  the  surface  (z  -  47 
m)  as  a  function  of  the  horizontal  distance  (solid  line).  Also  shown  are  the  sea 
surface  temperature  (dotted  line)  and  Initial  temperature  at  z  -  47  m  (dashed 
line).  The  temperature  drops  rapidly  as  expected  with  the  passage  of  the  gust 
front,  after  which  a  more  normal  value  is  reached.  Figure  7  indicates  cooling 


functions  of  height  at  (a)  x  =  217  km,  (b)  x  =  317  km,  and  (c)  x  =  417  km. 
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temperature  (•••).  initial  first  model  layer  temperature  ( — ),  and  predicted 
layer  temperature  ( — )  as  functions  of  the  horizontal  distance. 
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at  the  lowest  level  with  respect  to  the  initial  field.  Large  differences  between 
the  sea  surface  temperature  (T,  -  299.9  K)  and  the  first  model  layer  temperature 
will  significantly  affect  the  radiative  fluxes  at  lower  levels. 

3.3  Radiation  Budget  Diagnosis 

The  cloud  thermodynamic  and  microphysical  fields  simulated  by  the  CEM  are 
used  by  the  radiation  scheme  to  diagnose  the  radiative  budget  of  the  squall  line 
system.  For  diagnosis  of  solar  radiation  transfer,  a  solar  constant  of  1365  W 
m*^,  a  surface  albedo  of  0.05,  and  a  solar  zenith  angle  of  60 *  are  used. 

The  model  top  for  radiation  calculations  is  set  to  be  60  km  by  adding  two 
levels  above  the  CEM  domain.  These  two  levels  are  located  at  21  and  60  km. 
Temperature,  water  vapor,  and  ozone  profiles  at  21  km  are  assumed  to  be  the  same 
as  those  of  the  standard  tropical  atmosphere.  The  mixing  ratios  of  H2O  and  O3 
at  60  km  are  determined  by  equating  the  path  lengths  between  21  and  60  km  and 
those  obtained  from  the  detailed  tropical  profile.  The  temperature  at  60  km  is 
set  to  240  K,  which  is  the  average  temperature  between  21  and  60  km.  The 
differences  in  heating  rates  below  19  km  between  the  values  obtained  by  using  the 
detailed  profiles  above  the  cloud  model  domain  and  the  simplified  scheme  are 
within  0.02  K/day. 

The  average  IR  heating  rate  of  the  first  model  layer  (0-103.5  m)  is 
extremely  sensitive  to  the  surface  air  temperature  (T^) .  For  the  initial 
atmospheric  profiles,  the  IR  heating  rate  of  the  first  layer  is  -0.74  K/day, 
obtained  by  asstiming  that  the  surface  potential  temperature  and  water  vapor 
mixing  ratio  are  equal  to  those  at  47  m.  However,  the  heating  rate  is  -3.53 
K/day  if  the  surface  air  temperature  is  assumed  to  be  the  same  as  the  sea  surface 
temperature  and  if  the  surface  mixing  ratio  (qvo)  is  assumed  to  be  the  saturated 
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mixing  ratio.  Thus,  it  is  important  to  have  correct  Tj,  and  in  radiative 
calculations . 

Surface  fluxes  are  determined  from  the  flux-profile  relationships  developed 
by  Businger  et  al.  (1971)  and  Deardorff  (1972).  The  surface  potential 
temperature  flux,  (w'tf')o,  is  given  by 

-(w'  8’  )o  •  u.^,  »  -($,  -  6JC(z^)u,>  (3.1) 


where  S,  is  the  surface  value  of  S ,  -  ^(z„),  C  is  a  similarity  function  given 
below,  u.  is  the  friction  velocity,  and  z„  is  set  at  47  m  which  is  the  first 
model  level  height  above  the  surface.  The  similarity  function  for  the  unstable 
surface  layer  is 


where  y  -  (1  -  9$)^^*,  I  -  z/L,  k  is  the  Karman  constant  (0.35),  Zq  is  the 
roughness  length,  and  L  is  the  Monin-Obukhov  length.  From  Eq.  (3.1),  we  obtain 


K  ~ -iO,-  ejc(zj.  (3.3) 

A  detailed  profile  of  6  from  surface  to  z„  may  be  constructed  by  noting  that  6* 
is  independent  of  height.  Thus,  we  have 

«(z)  -  fl,  -  (^.  -  fl„)C(z„)/C(z). 

Similarly,  we  may  construct  a  detailed  profile  for  the  water  vapor  mixing  ratio, 
q^,  as  follows; 

q„(z)  -  q^,  -  (q,„  -  qy„)Ciz^) /C(z)  . 

In  the  present  case,  0,  is  the  potential  temperature  corresponding  to  a  constant 
sea  surface  temperature  of  299.9  K,  while  q^,  is  the  saturation  mixing  ratio 
corresponding  to  this  temperature;  6^  and  q^^  are  obtained  from  CEM.  From  GATE 
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data,  typical  values  for  and  L  are  0.6x10"^  m  and  -5  m,  respectively.  For  the 
initial  atmospheric  condition,  the  IR  heating  rate  for  the  first  model  layer  is 
-1.25  K/day  obtained  by  using  the  detailed  profiles  from  0  to  47  m  with  a 
vertical  resolution  of  1  m  in  the  radiation  calculation.  For  different 

(0.2xl0'^-0.6xl0'^  m)  and  L  (-1 — 20  m) ,  the  IR  heating  rate  ranges  from  -1.2  to 
-1.3  K/day. 

Based  on  numerical  experimentation,  by  setting 


Bo  -  B(z'). 

^vO  “  )  . 


(3.6) 


the  IR  heating  rates  computed  'r  m  using  the  GEM  model  resolution  agree  with 
those  from  the  detailed  ca’  alation  within  0.1  K/day.  The  parameter  z'  is 
related  to  L  by 


z'  -  0.2766  ♦  0.2324inlL|  .  (3.7) 

For  L  ranging  from  -1  to  -20,  z'  ranges  from  0.28  to  0.97  m. 


3.3.1  Heating  Rate  Fields 

In  the  following  radiation  calculations,  we  set  -  0.6x10'^  m,  and  set 
L  -  -5  m.  Figure  8  shows  the  detailed  cross-section  (x-z)  of  the  computed  solar 
heating  rates.  Solar  radiation  heats  the  upper  portions  of  the  clouds.  The 
maximum  heating  rate  is  -16.7  K/day,  which  occurs  at  the  top  of  the  deep 
cumulonimbi.  The  solar  heating  penetrates  through  the  anvil  top  for  -3  km.  The 
solar  heating  rate  field  in  a  squall  line  system  shows  a  significant  vertical  and 
horizontal  variability. 

Figure  9  shows  the  infrared  heating  rate  field  in  the  squall  line  system. 
Similar  to  solar  heating,  a  significant  spatial  variability  is  evident.  Strong 
infrared  cooling  occurs  at  the  cloud  top.  The  maximum  cooling  is  -24.9  K/day. 


a  contour  interval 


Cross  section  (x-z)  of  the  IR  heating  rate  with  a  contour  interval  of  2  K/day.  The  dashed 
contours  represent  cooling. 
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Little  infrared  heating  exists  near  the  cloud  base.  The  noticeable  heating  near 
the  surface  is  primarily  due  to  the  warmer  sea  surface  temperature ,  as  shown  in 
Fig.  7,  and  also  due  in  part  to  strong  absorption  of  rain.  The  maximum  infrared 
heating  is  -15.9  K/day,  which  is  smaller  than  the  maximum  cooling  value.  Like 
the  solar  heating  rate  field,  the  infrared  heating  rate  pattern  is  closely 
related  to  the  hydrometeor  mixing  ratio  field.  Comparing  Fig.  8  with  Fig.  9,  we 
see  that  the  maximum  solar  heating  occurs  below  the  maximum  infrared  cooling. 

Shown  in  Fig.  10  is  the  net  radiative  heating  rate.  Since  solar  heating 
and  IR  cooling  rates  are  comparable  in  the  upper  portion  of  the  anvil,  and  since 
the  maximum  solar  heating  occurs  below  the  maximum  infrared  cooling,  the  net 
heating  rate  shows  a  slight  cooling  at  the  cloud  top  and  a  slight  heating 
immediately  below.  At  the  cloud  top  where  water  phase  dominates,  IR  cooling  is 
much  higher  than  solar  heating,  resulting  in  strong  net  cooling.  Generally,  the 
net  radiative  heating  rate  field  shows  cooling  at  the  cloud  top  and  heating  at 
the  cloud  base  and  inside  clouds.  The  maximum  heating  and  cooling  rates  are  16.0 
and  18.0  K/day,  respectively.  Radiative  cooling  at  the  cloud  top  and  heating 
near  the  cloud  base  would  tend  to  generate  convective  mixing  of  the  cloud  layer. 

The  radiative  heating  rate  profiles  at  217,  317,  and  417  km  are  shown  in 
Fig.  11(a),  (b) ,  and  (c) ,  respectively.  In  these  cases,  since  cumulus/ 
cumulonimbus  clouds  are  connected  to  the  surface  through  falling  rain,  a  strong 
IR  heating  occurs  at  the  surface.  The  net  radiative  heating  profiles  at  all 
three  positions  exhibit  cooling  near  the  upper  portion  and  heating  within  the 
deeper  levels.  From  these  results,  it  is  clear  that  both  solar  and  IR  radiation 
are  important  in  determining  the  cloud  radiative  budget. 


3.3.2  Cloud  Radiative  Forcing 


section  (x-z)  of  the  net  (IR  plus  solar)  heating  rate  with  a  contour  interval  of  2  K/day 
ashed  contours  represent  cooling. 
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Fig.  11.  Solar  ( — ),  IR  (•  ■),  and  net  ( — )  heating  rates  as  functions  of  height  at  (a)  x  =  217  km,  (b)  x  = 
317  km,  and  (c)  x  =  417  km. 
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For  the  earth -atmosphere  system,  the  cloud  solar  radiative  forcing  is 
defined  as  the  difference  between  the  reflected  solar  fluxes  at  the  top  of  the 
atmosphere  for  clear  and  cloudy  skies.  The  cloud  IR  radiative  forcing  is  defined 
as  the  difference  between  the  outgoing  IR  fluxes  for  clear  and  cloudy  skies.  The 
net  cloud  radiative  forcing  is  the  sum  of  the  two.  The  concept  of  cloud 
radiative  forcing  provides  a  means  to  indicate  energy  gain/loss  of  the  earth- 
atmosphere  system  due  to  the  presence  of  clouds.  Cloud  radiative  forcing  may  be 
defined  with  reference  to  the  surface  and  the  atmosphere. 

Figure  12  shows  the  cloud  radiative  forcing  (for  the  earth -atmosphere 
system)  due  to  the  presence  of  a  squall  line.  As  noted  previously,  the  cloud  IR 
forcing  at  the  cloud  top  is  always  positive,  corresponding  to  the  heating  of  the 
system  due  to  the  cloud  greenhouse  effect,  while  the  cloud  solar  forcing  is 
always  negative,  indicating  the  cooling  of  the  system  by  the  cloud  albedo  effect. 
The  cloud  IR  forcing  ranges  from  14.7  to  162.4  W/m^,  while  the  solar  counterpart 
ranges  from  -129.3  to  -418.7  W/m^.  The  net  radiative  effect  of  clouds  leads  to 
the  cooling  of  the  earth -atmosphere  system.  It  is  Interesting  to  note  that  the 
maximum  IR  forcing  is  located  in  the  areas  of  deep  convection,  with  high  cloud 
tops  and  large  optical  depth.  The  net  cloud  radiative  forcing  at  the  top  of  the 
atmosphere  is  between  -109.0  and  -360.7  W/m^  in  a  squall  line  system. 

Figure  13  shows  the  cloud  radiative  forcing  at  the  surface.  The  IR  forcing 
is  always  positive  (23.4-50.2  W/m^)  because  the  cloud  base  emission  is  higher 
than  the  emission  from  a  clear  atmosphere  at  the  cloud  base  height.  From  Fig. 
13  it  is  noted  that  the  negative  surface  cloud  forcing  values  due  to  solar 
radiation  (-134.4 — 443.8  W/m^)  are  substantially  similar  to  that  presented  in 
Fig.  12  (dashed  line).  This  is  because  the  atmosphere  is  largely  transparent 
with  respect  to  solar  radiation.  The  net  surface  cloud  forcing  ranges  between 
-111.0  and  -397.9  W/m^- 


mgs  at  the  top  ot  the  atmosphere  tor  solar 
functions  of  the  horizontal  distance. 


Surface 
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Figure  14  depicts  the  cloud  radiative  forcing  for  the  atmosphere,  which  is 
the  difference  between  the  cloud  radiative  forcing  at  the  top  and  that  at  the 
surface.  The  solar  warming  due  to  the  presence  of  clouds  is  small,  ranging  from 
3.4  to  30.9  W/m^.  The  cloud  IR  forcing  is  also  positive  in  the  areas  where  high 
clouds  are  present:  a  large  value  of  112.2  W/m^  in  the  case  of  deep  cumulonimbi 
is  shown.  For  low  clouds,  the  cloud  IR  forcing  is  negative  with  values  up  to 
-18.9  W/m^.  The  net  atmospheric  cloud  radiative  forcing  is  between  -3.1  and 
138.9  W/m*  in  the  present  study. 

Referring  to  Figs.  12,  13,  and  14,  we  conclude  that  in  a  squall  line  system 
(or  tropical  deep  convective  areas) ,  the  cloud  solar  forcing  for  the  earth- 
atmosphere  system  is  largely  confined  to  the  sea  surface,  while  the  cloud  IR 
forcing  is  mainly  confined  within  the  atmosphere.  The  net  radiation  effect  is 
cooling  at  the  surface  and  heating  within  the  atmosphere.  Thus,  the  radiative 
effect  of  clouds  is  similar  to  the  latent  heat  release  in  that  both  remove  heat 
from  the  surface  and  deposit  it  into  the  atmosphere  (Ramanthan,  1987).  However, 
the  surface  cooling  resulting  from  the  reduction  of  solar  radiation  is 
significantly  larger  than  the  heating  in  the  atmosphere  produced  by  IR  radiative 
exchanges.  The  presence  of  clouds  significantly  modifies  the  vertical 
distribution  of  radiative  heating/cooling  which  would  significantly  affect 
atmospheric  convection  and  dynamic  processes. 


Section  4 


SUMMARY 

We  have  interfaced  the  newly  developed  radiation  scheme  with  a  CEM  to 
Investigate  the  interactions  between  the  radiation  field  and  tropical  convective 
systems.  The  thermodynamic,  water  vapor,  and  bulk  hydrometeor  fields  provided 
by  the  CEM  are  used  as  diagnostic  inputs  for  radiation  calculations.  The 
radiative  properties  of  hydrometeors,  including  water  droplets,  ice  crystals, 
rain,  snow,  and  graupel  have  been  treated  explicitly.  The  single-scattering 
properties  of  ice  crystals  are  parameterized  as  functions  of  mean  effective  size 
and  ice  water  content,  based  on  a  light  scattering  program  for  hexagonal  plates 
and  colvimns.  The  scattering  and  absorption  properties  of  water  droplets  are 
represented  as  functions  of  mean  effective  radius  and  liquid  water  content  based 
on  Mie  scattering  calculations.  Moreover,  the  radiative  properties  of  rain, 
snow,  and  graupel  are  calculated  from  Mie  theory  using  Marshall -Palmer 
distributions  that  are  employed  in  the  bulk  microphysical  parameterization.  A 
(-four-stream  radiative  transfer  scheme  has  been  developed  for  flux  calculations 
in  both  solar  and  Infrared  spectra.  For  nongray  gaseous  absorption  due  to  HjO, 
CO2,  O3,  CH*,  and  N2O,  the  correlated  k-distribution  method  is  used,  which  can  be 
readily  Incorporated  into  scattering  models. 

The  initial  thermodynamic  state  based  on  the  GATE  Phase -111  mean  sounding 
is  used  along  with  the  present  CEH  to  simulate  a  tropical  squall  line  system. 
The  CEM  was  able  to  simulate  many  characteristic  features  of  tropical  squall 
lines,  such  as  the  size  and  shape  of  the  clouds  associated  with  the  tropical 
mesoscale  convective  system.  The  radiation  budgets  are  calculated  by  using  the 
CEM  simulated  thermodynamic  and  hydrometeor  fields. 

By  examining  the  radiative  heating  rate  fields,  it  is  found  that  solar 
radiation  heats  the  upper  part  of  the  cloud,  which  would  stabilize  the  cloud 
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layer  and  evaporate  the  cloud  particles.  On  the  contrary,  IR  cooling  occurs  at 
the  cloud  top,  which  would  enhance  the  convective  instability.  Significant 
heating  at  the  cloud  base  does  not  take  place  in  deep  anvils  because  of  the 
presence  of  low  clouds  and  because  of  the  high  cloud-base  temperature  and  large 
air  density.  Because  the  maximuim  solar  heating  occurs  below  the  maximum  IR 
cooling,  the  net  radiative  heating  shows  a  pattern  of  cooling  above  heating  in 
the  upper  part  of  the  cloud  layer.  This  pattern  could  result  in  convection  at 
the  cloud  top  and  promote  entrainment. 

The  tropical  mesoscale  system  has  a  significant  Impact  on  the  radiative 
budget  of  the  earth  atmosphere  system.  The  system  reduces  the  loss  of  IR  fluxes 
emitted  from  the  top  of  the  atmosphere  and  increases  IR  emission  to  the  surface. 
In  the  deep  convective  areas,  the  reduction  of  the  loss  of  IR  fluxes  emitted  from 
the  top  of  the  atmosphere  can  be  as  large  as  160  W/m^,  which  is  redistributed  to 
the  surface  and  the  atmosphere  by  -50  and  -110  W/m^,  respectively.  The  tropical 
mesoscale  system  significantly  Increases  the  solar  albedo.  The  deep  convective 
clouds  reflect  -420  W/m*  more  flux  than  clear  sky  and  largely  reduce  the  solar 
fluxes  available  at  the  sea  surface. 

The  radiative  fields,  including  both  heating  rate  and  cloud  radiative 
forcing,  show  significant  vertical  and  horizontal  variabilities  in  the  presence 
of  the  tropical  mesoscale  systems.  The  circulation  patterns  can  be  profoundly 
modulated  by  these  variabilities.  It  should  be  noted  that  the  vertical  grid  used 
in  the  present  study  is  stretched  from  -100  m  near  the  surface  to  -1000  m  at  the 
top  of  the  CEM.  Thus,  the  vertical  resolution  at  the  deep  convective  cloud  top 
for  radiation  calculations  is  -1  km.  This  is  not  sufficient  to  resolve  the 
small-scale  vertical  structure  involving  the  heating  rate  profile  at  the  cloud 
top  for  the  investigation  of  the  interactions  of  radiation,  microphysical 


processes,  and  dynamic  motions.  The  effects  of  vertical  resolution  on  these 
interactions  require  further  numerical  studies. 
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ABSTRACT 

A  new  approach  for  parameterization  of  the  broadband  solar  and  infrared 
radiative  properties  of  ice  clouds  has  been  developed.  This  parameterization 
scheme  integrates  in  a  coherent  manner  the  5 -four-stream  approximation  for 
radiative  transfer,  the  correlated  k-distribution  method  for  nongray  gaseous 
absorption,  and  the  scattering  and  absorption  properties  of  hexagonal  ice 
crystals.  We  use  a  mean  effective  size,  representing  an  area  weighted  mean 
crystal  width,  to  account  for  the  ice  crystal  size  distribution  with  respect  to 
radiative  calculations.  Based  on  physical  principles,  the  basic  single¬ 
scattering  properties  of  ice  crystals,  including  the  extinction  coefficient 
divided  by  ice  water  content,  single -scattering  albedo,  and  expansion 
coefficients  of  the  phase  function,  can  be  parameterized  using  third-degree 
polynomials  in  terms  of  the  mean  effective  size.  In  the  development  of  this 
parameterization,  we  use  the  results  computed  from  a  light  scattering  program 
that  includes  a  geometric  ray-tracing  program  for  size  parameters  larger  than  30 
and  the  exact  spheroid  solution  for  size  parameters  less  than  30.  The 
computations  are  carried  out  for  11  observed  ice  crystal  size  distributions  and 
cover  the  entire  solar  and  thermal  infrared  spectra.  Parameterization  of  the 
single-scattering  properties  is  shown  to  provide  an  accuracy  within  about  1%. 
Comparisons  have  been  carried  out  between  results  computed  from  the  model  and 
those  obtained  during  the  1986  cirrus  FIRE  IFO.  We  show  that  the  model  results 
can  be  used  to  reasonably  interpret  the  observed  IR  emissivities  and  solar  albedo 
involving  cirrus  clouds.  The  newly  developed  scheme  has  been  employed  to 
investigate  the  radiative  effects  of  ice  crystal  size  distributions.  For  a  given 
ice  water  path,  cirrus  clouds  with  smaller  mean  effective  sizes  reflect  more 
solar  radiation,  trap  more  infrared  radiation,  and  produce  stronger  cloud-top 
cooling  and  cloud-base  heating.  The  latter  effect  would  enhance  the  in-cloud 


ii 


heating  rate  gradients.  Further,  the  effects  of  ice  crystal  size  distribution 
in  the  context  of  IR  greenhouse  versus  solar  albedo  effects  involving  cirrus 
clouds  are  presented  with  the  aid  of  the  upward  flux  at  the  top  of  the 
atmosphere.  In  most  cirrus  cases,  the  IR  greenhouse  effect  outweighs  the  solar 
albedo  effect.  One  exception  occurs  when  a  significant  number  of  small  ice 
crystals  is  present.  The  present  scheme  for  radiative  transfer  in  the  atmosphere 
involving  cirrus  clouds  is  well  suited  for  incorporation  in  numerical  models  to 
study  the  climatic  effects  of  cirrus  clouds,  as  well  as  to  investigate 
interactions  and  feedbacks  between  cloud  microphysics  and  radiation. 
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1 .  Introduction 

Cirrus  clouds  are  globally  distributed,  being  present  at  all  latitudes  and 
without  respect  to  land  or  sea  or  season  of  the  year.  They  regularly  cover  about 
20-30%  of  the  globe  and  strongly  influence  weather  and  climate  processes  through 
their  effects  on  the  radiation  budget  of  the  earth  and  the  atmosphere  (Liou, 
1986)  .  The  importance  of  cirrus  clouds  in  weather  and  climate  research  can  be 
recognized  by  the  intensive  field  observations  that  have  been  conducted  as  a 
major  component  of  the  First  ISCCP  Regional  Experiment  in  October-November  1986 
(Starr,  1987)  and  more  recently  in  November  1991. 

Cirrus  clouds  possess  a  number  of  unique  features.  In  addition  to  being 
global  and  located  high  in  the  troposphere  and  extending  to  the  lower 
stratosphere  on  some  occasions,  they  contain  almost  exclusively  nonspherical  ice 
crystals  of  various  shapes ,  such  as  bullet  rosettes ,  plates ,  and  colvimns .  There 
are  significant  computational  and  observational  difficulties  in  determining  the 
radiative  properties  of  cirrus  clouds.  A  reliable  and  efficient  determination 
of  the  radiative  properties  of  cirrus  clouds  requires  the  fundamental  scattering 
and  absorption  data  involving  nonspherical  ice  crystals.  In  addition, 
appropriate  incorporation  of  gaseous  absorption  in  scattering  cloudy  atmospheres 
and  an  efficient  radiative  transfer  methodology  are  also  required. 

Although  parameterization  of  the  broadband  radiative  properties  for  cirrus 


clouds  has  been  presented  by  Liou  and  Wlttman  (1979)  in  terms  of  ice  water 
content,  such  parameterization  used  the  scattering  and  absorption  properties  of 
circular  cylinders  without  accounting  for  the  effects  of  the  hexagonal  structure 
of  ice  crystals.  Moreover,  the  effects  of  ice  crystal  size  distribution  were  not 
included  in  the  parameterization.  Using  area-equivalent  or  volume -equivalent  ice 
spheres  to  approximate  hexagonal  ice  crystals  for  scattering  and  absorption 
properties  has  been  shown  to  be  inadequate,  and  frequently  misleading  This  is 
evident  in  interpreting  the  scattering  and  polarization  patterns  from  ice  clouds 


2 


(Takano  and  Liou,  1989)  and  the  observed  radiative  properties  of  cirrus  clouds, 
in  particular,  the  cloud  albedo  (Stackhouse  and  Stephens,  1991). 

In  this  paper,  we  wish  to  develop  a  new  approach  for  the  parameterization 
of  the  broadband  solar  and  infrared  radiative  properties  of  ice  crystal  clouds. 
Three  major  components  are  integrated  in  this  parameterization,  including  the 
scattering  and  absorption  properties  of  hexagonal  ice  crystals,  the  ^-four-stream 
approximation  for  radiative  transfer,  and  the  correlated  k-distribution  method 
for  nongray  gaseous  absorption.  Compared  with  more  sophisticated  models  and 
aircraft  observations ,  the  present  parameterization  has  been  shown  to  be  accurate 
and  efficient  for  flux  and  heating  rate  calculations.  In  Section  2,  we  present 
parameterization  of  the  single -scattering  properties  of  ice  crystals.  The  manner 
in  which  the  radiative  flux  transfer  is  parameterized  using  the  i -four-stream 
approximation  and  the  correlated  k-distribution  is  discussed  in  Section  3. 
Section  4  presents  some  comparisons  between  theoretical  results  and  observed  data 
involving  cloud  emissivity  and  albedo.  In  Section  5,  we  present  the  effect  of 
ice  crystal  size  distribution  on  cloud  heating  rates  and  on  the  question  of  cloud 
radiative  forcing.  Summary  and  conclusions  are  given  in  Section  6, 

2.  Parameterization  of  the  Single-Scattering  Properties  of  Ice  Crystals 
a.  Physical  Bases 

The  calculations  of  the  single -scattering  properties,  including  the  phase 
function,  single -scattering  albedo  and  extinction  coefficient,  require  a  light 
scattering  program  and  the  detailed  particle  size  distribution.  The  calculations 
are  usually  time  consuming.  If  radiation  calculations  are  to  interact  with  an 
evolving  cloud  where  particle  size  distribution  varies  as  a  function  of  time 
and/or  space,  the  computer  time  needed  for  examining  just  this  aspect  of  the 
radiation  program  would  be  formidable,  even  with  a  super  computer.  Thus,  there 
is  a  practical  need  to  simplify  the  computational  procedure  for  the  calculation 


of  the  single-scattering  properties  of  cloud  particles.  Since  spheres  scatter 
an  amount  of  light  proportionate  to  their  cross-section  area,  a  mean  effective 
radius,  which  is  defined  as  the  mean  radius  that  is  weighted  by  the  cross-section 
area  of  spheres,  has  been  used  in  conjunction  with  radiation  calculations  (Hansen 
and  Travis,  1974).  Higher  order  definition,  such  as  dispersion  of  the  droplet 
size,  may  be  required  in  order  to  more  accurately  represent  the  droplet  size 
distribution. 

Ice  crystals  are  nonspherical  and  ice  crystal  size  distributions  are 
usually  expressed  in  terms  of  the  maximum  dimension  (or  length)  .  Representation 
of  the  size  distribution  for  ice  crystals  is  much  more  involved  than  that  for 
water  droplets.  To  the  extent  that  scattering  of  light  is  proportional  to  the 
cross-section  area  of  nonspherical  particles,  we  may  use  a  mean  effective  size 
analogous  to  the  mean  effective  radius  defined  for  spherical  water  droplets  as 
follows: 

D,  -  r^“D  •  DLn(L)dL/  r^"'DLn(L)  dL,  (2.1) 

J^iin  J^ln 

where  D  is  the  width  of  an  ice  crystal,  n(L)  denotes  the  ice  crystal  size 
distribution,  and  and  are  the  minimtam  and  maximum  lengths  of  ice 
crystals,  respectively.  A  similar  definition  for  circular  cylinders  has  also 
been  proposed  by  Platt  and  Harshvardhan  (1988).  Based  on  aircraft  observations 
by  Ono  (1969)  and  Auer  and  Veal  (1970),  the  width  may  be  related  to  the  length 
L.  It  follows  that  the  mean  effective  width  (or  size)  can  be  defined  solely  in 
terms  of  ice  crystal  size  distribution.  The  geometric  cross  section  area  for 
oriented  hexagonal  ice  crystals  generally  deviates  from  DL  [see  Eq.  (2.4)  for  the 
condition  of  random  orientation]  .  To  the  extent  that  D  is  related  to  L,  the 
definition  of  D,  in  Eq.  (2.1),  which  is  an  approach  to  represent  ice  crystal  size 
distribution,  should  be  applicable  to  ice  crystals  with  hexagonal  structure.  The 
nximerator  in  Eq.  (2.1)  is  related  to  the  ice  water  content  (IWC)  in  the  form 
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IWC  =  Pi  [*"“0  •  DLn(L)dL, 

O  J^nin 

where  the  volume  of  a  hexagonal  ice  crystal,  3  JZ  D^L/8,  is  used  and  is  the 
density  of  ice.  We  shall  confine  our  study  to  the  use  of  the  mean  effective  size 
to  represent  ice  crystal  size  distribution  in  the  single -scattering  calculation 
for  ice  crystals . 

The  extinction  coefficient  is  defined  by 


^  *  r^"a(D,L)n(L)  dL.  (2.3) 
JLmln 

where  a  is  the  extinction  cross  section  for  a  single  crystal.  In  the  limits  of 
geometric  optics  and  using  hexagonal  ice  crystals  that  are  randomly  oriented  in 
space,  the  extinction  cross  section  may  be  expressed  by  (Takano  and  Liou,  1989) 


a 


D  +  Lj 


(2.4) 


Substituting  Eq.  (2.4)  into  Eq.  (2.3)  and  using  the  definitions  of  D,  and  IWC, 
we  have 


IWC  • 


fW  D2n(L)  jjL/  D2Ln(L)  dL  +  -A_  ^ 
Pi  Kin  Kin  v/3pi  ^ 


(2.5a) 


The  first  term  on  the  right  cannot  be  defined  in  terms  of  D,  directly.  However, 
since  D  <  L,  this  term  should  be  much  smaller  than  the  second  term,  which  also 
in'/^olves  a  factor  4/73.  To  the  extent  that  D  is  related  to  L,  the  first  term  may 
be  approximated  by  a  +  b'/D, ,  where  b'  «  4/73pi  and  a  is  a  certain  constant. 
Using  this  argument,  we  obtain  the  following  relationship: 


/9  ^  IWC(a  +  b/D,), 


(2.5b) 


where  b  -  b'  +  Using  the  observed  ice  crystal  size  distributions 

described  below,  we  also  show  that  ^/IWC  and  1/D,  are  linearly  related  in  the 
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solar  spectral  region.  Based  on  the  preceding  analysis,  it  is  clear  that  the 
extinction  coefficient  is  a  function  of  both  IWC  and  mean  effective  size. 

Because  cloud  absorption  is  critically  dependent  on  the  variation  of  the 
single-scattering  albedo,  it  must  be  accurately  parameterized.  For  a  given  ice 
crystal  size  distribution,  the  single -scattering  albedo,  w,  is  defined  by 

1  -  w  -  r^a.n(L)  dL/ f^an(L)  dL,  (2.6) 

J^Bln  J^aln 

where  o,  denotes  the  absorption  cross  section  for  a  single  crystal.  When 
absorption  is  small,  o,  is  approximately  equal  to  the  product  of  the  imaginary 
part  of  the  refractive  index  of  ice,  m^,  and  the  particle  volume,  viz., 

^  _  3^3  xm^  (A)  ^2^  (2.7) 

a.  - ^ - D  L, 

where  A  is  the  wavelength.  Using  the  extinction  and  absorption  cross  sections 
defined  in  Eqs.  (2.4)  and  (2.7)  and  noting  that  D  is  related  to  L,  based  on 
observations,  we  obtain 

1  -  w  »  c  +  dD, ,  (2.8) 

where  c  and  d  are  certain  coefficients. 

In  the  preceding  discussion,  we  have  used  the  geometric  optics  limit  to 
derive  the  expressions  for  the  extinction  coefficient  and  single -scattering 
albedo.  In  view  of  the  observed  ice  crystal  sizes  in  cirrus  clouds  (-20-2000 
nm)  ,  the  simple  linear  relationships  denoted  in  Eqs.  (2.5b)  and  (2.8)  should  be 
valid  for  solar  wavelengths  (0.2-4  tm)  .  For  thermal  infrared  wavelengths  (e.g.  , 
10  pm) ,  the  geometric  optics  approximation  may  not  be  appropriate  for  small  ice 
crystals.  However,  we  note  from  aircraft  observations  that  there  is  a  good 
linear  relationship  between  the  extinction  coefficient  in  the  infrared  spectrum 
and  the  extinction  coefficient  derived  based  upon  the  large -particle 
approximation  (Foot,  1988). 
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b.  A  Generalized  Parameterization 

The  linear  relationship  between  fi/lMC  and  1/D,  shown  in  Eq.  (2.5b)  is 
derived  based  on  the  geometric  optics  approximation  and  the  assumption  that  ice 
crystals  are  randomly  oriented  in  space.  The  linear  relationship  between  w  and 
D,  shown  in  Eq.  (2.8)  is  based  on  the  asstimption  that  ice  crystal  absorption  is 
small  and  that  ice  crystals  are  randomly  oriented.  For  general  cases,  we  would 
expect  that  higher-order  expansions  may  be  needed  to  define  more  precisely  the 
single -scattering  properties  of  ice  crystals  in  terms  of  the  mean  effective  size. 
Thus  we  postulate  that 


0  «■  IWC 


(2.9) 

(2.10) 


where  a^  and  b„  are  certain  coefficients,  which  must  be  determined  from  numerical 
fitting,  and  N  is  the  total  number  of  terms  required  to  achieve  a  prescribed 
accuracy.  When  N«l,  Eqs.  (2.9)  and  (2.10)  are  exactly  the  same  as  Eqs.  (2.5b) 
and  (2.8).  When  N-2,  the  term,  1/DJ  is  proportional  to  the  variance  of  ice 
crystal  size  distribution.  Based  on  numerical  experimentation  described  in 
subsection  c,  we  find  N-2  is  sufficient  for  the  extinction  coefficient  expression 
to  achieve  an  accuracy  within  1%.  For  the  single-scattering  albedo,  we  find  that 
N-3  is  required. 

For  nonspherical  particles  randomly  oriented  in  space,  the  phase  function 
is  a  function  of  the  scattering  angle,  0.  The  phase  function  is  usually  expanded 
in  a  series  of  Legendre  polynomials  P,  in  radiative  transfer  calculations  in  the 
form 

M 

P(cos  0)  »  (cos  0)  ,  (2.11) 


where  we  set  wq  ”  1.  Since  the  phase  function  is  dependent  on  ice  crystal  size 
distribution,  the  expansion  coefficients  must  also  be  related  to  ice  crystal  size 
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distribution,  which  is  represented  by  the  mean  effective  size  in  the  present 
study.  In  the  case  of  hexagonal  ice  crystals,  in  addition  to  the  diffraction, 
scattered  energy  is  also  produced  by  the  5-function  transmission  through  parallel 
planes  at  0  -  0  (Takano  and  Liou,  1989).  Using  the  similarity  principle  for 
radiative  transfer,  the  expansion  coefficients  in  the  context  of  four-stream 
approximation  can  be  expressed  Ly 

w,  =  (1  -  fi)«;  +  (2i  ♦  1),  -e  =  l,2,3,4,  (2.12) 


where  (j,*  represents  the  expansion  coefficients  for  the  phase  function  in  which 
the  forward  6 -function  peak  has  been  removed,  and  f^  is  the  contribution  from  the 
forward  5-function  peak.  In  our  notation,  -  3g,  where  g  is  the  asymmetry 
factor.  The  5 -function  peak  contribution  has  been  evaluated  by  Takano  and  Liou 
and  is  a  function  of  ice  crystal  size.  The  f^  value  increases  with  increasing 
L/D  value  due  to  a  greater  probability  for  plane -parallel  transmission.  Ue  may 
express  wj  and  f^  in  terms  of  the  mean  effective  size  as  follows; 


(2.13a) 


E 

n*0 


(2.13b) 


where  c„  |  and  d„  are  certain  coefficients.  Based  on  numerical  experimentation 
described  in  subsection  c,  we  find  N-3  is  sufficient  to  achieve  an  accuracy 
within  1%. 

In  the  thermal  infrared  wavelengths,  halo  and  5 -transmission  peak  features 
in  the  phase  function  are  largely  suppressed  due  to  strong  absorption.  For  this 
reason  and  to  a  good  approximation,  we  may  use  the  asymmetry  factor  to  represent 
the  phase  function  via  the  Henyey-Greenstein  function  in  the  form 
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w,  =  (2i  +1)  g*.  (2.14) 

The  asynunetry  factor  may  also  be  expressed  In  terms  of  the  effective  mean  size 
as  follows: 

g(IR)  =  £  c^  D.",  <2.15) 

n«0 

where  c^  again  is  a  certain  coefficient  and  N-3  is  sufficient  in  the  expansion, 
c.  Determination  of  the  Coefficients  in  the  Parameterization 

The  coefficients  in  Eqs.  (2.9),  (2.10),  (2.13),  and  (2.15)  are  determined 
from  numerical  fitting  to  the  data  computed  from  "exact”  light  scattering  and 
absorption  programs  that  include  the  following.  For  size  parameters  larger  than 
about  30,  we  use  the  results  computed  from  a  geometric  ray- tracing  program  for 
hexagonal  ice  crystals  developed  by  Takano  and  Liou  (1989) .  Note  that  the 
parameterization  coefficients  may  be  updated  if  the  results  for  other  t3T>es  of 
ice  crystals,  such  as  bullet  rosettes  and  hollow  columns,  are  available.  For 
size  parameters  less  than  about  30,  the  laws  of  geometric  optics  are  generally 
not  applicable.  Since  the  exact  solution  for  hexagonal  ice  crystals  based  on 
either  the  wave  equation  approach  or  other  integral  methods  for  size  parameters 
on  the  order  of  20-30  has  not  been  developed,  we  use  the  results  computed  from 
a  light  scattering  program  for  spheroids  developed  by  Asano  and  Sato  (1980) .  We 
have  tested  and  modified  this  program  (Takano  et  al.,  1992)  and  found  that  it  can 
be  applied  to  size  parameters  less  than  about  30.  Naturally  occurring  ice 
crystals  will  have  a  hexagonal  structure  that  cannot  be  approximated  by 
spheroids.  However,  small  size  parameters  are  usually  associated  with  infrared 
wavelengths,  where  ice  is  highly  absorbing.  It  is  likely  that  the  detailed  shape 
factor  may  not  be  critical  in  scattering  and  absorption  calculations.  Further 
inquiry  into  this  problem  appears  necessary. 

In  the  present  "exact"  computations,  11  ice  crystal  size  distributions  from 
in-situ  aircraft  observations  were  employed.  Table  1  shows  the  number  densities, 
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IWCs,  and  mean  effective  sizes  for  these  11  size  distributions.  The  first  two 
t3T>es  are  for  the  cirrostratus  (Cs)  and  cirrus  uncinus  presented  by  He3nnsfield 
(1975),  while  the  next  two  types  are  modified  distributions  given  by  Heymsfield 
and  Platt  (1984)  corresponding  to  warm  and  cold  cirrus  clouds.  The  size 
distributions  of  cirrus  clouds  at  temperatures  of  -20,  -40,  and  -60”C  were 
obtained  from  the  parameterization  results  presented  by  Heymsfield  and  Platt 
(1984).  Other  cirrus  types,  Ci  (Oct.  22,  Oct.  25,  Nov.  1,  and  Nov.  2),  were  the 
ice  crystal  size  distributions  derived  from  the  1986  FIRE  cirrus  experiments 
(Heymsfield,  personal  communication).  The  11  size  distributions  cover  a 
reasonable  range  of  cloud  microphysical  properties  in  terms  of  IWC  (6.6  x  10“*  - 
0.11  gm“®)  and  mean  effective  size  (23.9  -  123.6  pm)  as  shown  in  Table  1.  For 
scattering  and  absorption  calculations,  these  size  distributions  have  been 
discretized  in  five  regions.  The  aspect  ratios,  L/D,  used  are  20/20,  50/40, 
120/60,  300/100,  and  750/160  in  units  of  pm/pm,  roughly  corresponding  to  the 
observations  reported  by  Ono  (1969)  and  Auer  and  Veal  (1970). 

In  the  single-scattering  calculations,  the  refractive  indices  for  ice, 
compiled  by  Warren  (1984) ,  were  used.  In  order  to  resolve  the  variation  in  the 
refractive  index  of  ice  and  account  for  the  gaseous  absorption,  six  and  12  bands 
were  selected  for  solar  and  thermal  IR  regions,  respectively.  The  spectral 
division  is  shown  in  Table  2,  where  the  index  for  the  spectral  band  (i  - 
1,2. . . ,18)  is  defined.  The  complex  indices  of  refraction  at  each  wavelength  are 
averaged  values  over  the  spectral  band,  weighted  by  the  solar  Irradiance 
(Thekaekara,  1973)  in  solar  spectrum  and  by  the  Planck  function  (T  -  -40*C)  in 
thermal  IR  spectrum.  This  temperature  is  a  typical  value  for  cirrus  clouds.  The 
single-scattering  calculations  were  carried  out  at  the  central  wavelength  for 
each  spectral  band.  The  coefficients  in  Eqs.  (2.9),  (2.10),  (2.13),  and  (2.15) 
determined  by  numerical  fitting  using  the  "exact"  results  are  listed  in  Tables 
3-5.  In  the  solar  spectrum  where  the  law  of  geometric  optics  is  valid,  the 


10 


extinction  coefficients  for  a  given  size  distribution  are  the  same  regardless 
of  the  wavelength,  as  shown  in  Table  3.  The  optical  depth  r  -  ^Az  can  be 
obtained  if  the  cloud  thickness  Az  is  given.  Figures  1-2  show  the  fitting  for 
;9/IWC,  Z  and  g  in  the  spectral  intervals  1.9 -2. 5  urn  and  800-980  cm“^, 
respectively.  For  the  solar  spectrum,  we  use  band  number  four  (1.9 -2. 5  nm)  as 
an  example,  because  the  single -scattering  properties  in  this  band  significantly 
depend  on  size  distribution.  The  linear  relationship  between  j8/IWC  and  l/D, 
postulated  in  Eq.  (2.5b)  is  clearly  shown  in  Fig.  la.  This  linear  relationship 
is  valid  for  all  other  solar  bands ,  because  the  extinction  cross  sections  for 
solar  bands  in  the  limit  of  the  geometric  optics  approach  are  the  same,  as 
pointed  out  previously.  The  nonlinearity  between  u>  and  D,  becomes  important  for 
small  Z.  This  is  because  the  linear  relationship  developed  in  Eq.  (2.8)  is 
based  on  the  weak  absorption  assumption.  The  band  800-980  cm“^  is  located  in  the 
atmospheric  window,  where  the  greenhouse  effect  of  clouds  is  most  pronounced. 
It  can  be  seen  that  very  good  fits  for  single-scattering  properties  in  both  solar 
and  infrared  spectra  are  obtained.  The  relative  errors  are  less  than  -1%.  Other 
spectral  bands  also  show  similar  accuracies. 

It  is  known  that  the  minimum  ice  crystal  length  that  the  present  optical 
probe  can  measure  is  about  20  /im.  To  investigate  the  potential  effects  of  small 
ice  crystals  that  may  exist  in  cirrus  clouds  on  parameterizations  of  the  single - 
scattering  properties,  we  used  the  11  observed  ice  crystal  size  distributions  and 
extrapolated  these  distributions  from  20  to  10  ^m  in  the  logarithmic  scale.  The 
resulting  mean  effective  sizes  range  from  18.9  to  122.7  ftm.  Using  these  sizes, 
the  numerical  fitting  coefficients  vary  only  slightly  and  do  not  affect  the 
accuracy  of  the  preceding  parameterized  single-scattering  properties.  In  the 
present  study,  ice  crystals  are  assumed  to  be  randomly  oriented  in  space.  In 
cases  when  cirrus  clouds  contain  horizontally  oriented  ice  crystals, 
parameterizations  would  require  modifications.  The  extinction  coefficient. 
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single-scattering  albedo  and  phase  function  would  depend  on  the  incident 
direction.  The  radiative  transfer  scheme  would  also  require  adjustments  to 
account  for  anisotropic  properties.  Since  fluxes  are  involved  in  the 
parameterization,  we  anticipate  that  the  modifications  should  not  significantly 
affect  the  results  and  conclusions  derived  from  this  study.  At  any  rate,  this 
is  an  area  requiring  further  research  efforts. 

3.  Parameterization  of  Radiative  Flux  Transfer 
a.  Radiative  Transfer  Scheme 

We  follow  the  6-four-stream  model  developed  by  Liou  et  al.  (1988)  for  the 
calculations  of  solar  radiative  flux  transfer  in  a  single  homogeneous  layer.  The 
solution,  like  various  two-stream  methods,  is  in  analytic  form  so  that  the 
computational  effort  involved  is  minimal.  As  demonstrated  in  that  paper,  results 
from  the  6 -four- stream  approximation  can  yield  relative  accuracies  within  -5%. 

To  obtain  a  single  treatment  of  solar  and  infrared  radiation,  we  extend  the 
6 -four -stream  approach  to  the  transfer  of  infrared  radiation,  in  which  the 
optical  depth  dependence  of  the  Planck  function  must  be  known.  The  Planck 
function  can  be  approximately  expressed  in  terms  of  the  optical  depth  in  the  form 

B(t)  =  Boexp(r/riinBi/Bo)  ,  (3.1) 
where  is  the  optical  depth  of  a  layer,  and  Bq  and  Bj  are  the  Planck  functions 
corresponding  to  the  temperature  at  the  top  and  bottom  of  this  layer, 
respectively.  This  exponential  approximation  for  the  Planck  function  in  optical 
depth  has  an  accuracy  similar  to  the  linear  approximation  for  the  Planck  function 
developed  by  Wiscombe  (1976) .  Since  the  direct  solar  radiation  source  has  an 
exponential  function  form  in  terms  of  optical  depth,  the  formulation  of  the  6- 
four-stream  approximation  for  infrared  wavelengths  is  the  saime  as  that  for  solar 
wavelengths.  For  this  reason,  the  computer  program  is  simplified  and  the 
computational  speed  is  enhanced. 
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In  the  limit  of  pure  molecular  absorption  at  infrared  wavelengths,  Tj  may 
approach  zero  for  some  bands.  The  major  advantage  of  using  Eq.  (3.1)  is  that  the 
solution  of  the  transfer  equation  is  numerically  stable  when  -»  0.  If  the 
linear  approximation  in  r  is  used,  the  round-off  error  may  lead  to  an  unstable 
solution  when  is  very  small.  We  note  that  the  higher  order  polynomial 
approximation  in  r  for  the  Planck  function  can  also  lead  to  an  unstable  solution. 
As  in  the  case  of  solar  wavelengths,  we  perform  relative  accuracy  checks  for 
infrared  wavelengths  with  respect  to  "exact"  results  computed  from  the  adding 
method  for  radiative  transfer  (Liou,  1992).  Numerical  results  reveal  that  the 
errors  in  the  j-four-stream  scheme  in  computing  IR  cloud  emissivity  are  typically 
less  than  1%.  Systematic  equations  for  this  scheme  that  are  required  in 
numerical  calculations  are  given  in  the  Appendix. 

For  applications  of  the  S -four-stream  scheme  to  nonhomogeneous  atmospheres, 
the  atmosphere  is  divided  into  N  homogeneous  layers  with  respect  to  the  s ingle - 
scattering  albedo  and  phase  function.  The  4  x  N  unknown  coefficients  in  the 
analytic  solution  for  the  transfer  equation  are  determined  following  the 
procedure  described  in  Liou  (1975).  A  numerically  stable  program  has  been 
developed  to  solve  the  system  of  linear  equations  by  utilizing  the  property  that 
the  coefficient  matrix  is  a  sparse  matrix.  While  this  method  is  similar  to  the 
Gaussian  elimination  method  with  back  substitution,  it  is  carefully  optimized  so 
as  to  minimize  the  mathematical  operations.  We  perform  partial  pivoting 
(interchanging  rows  only)  to  guarantee  that  a  solution  always  exists.  In  order 
to  avoid  a  fatal  overflow  problem,  we  follow  the  scaling  technique  developed  by 
Stamnes  and  Conklin  (1984)  so  that  numerical  solutions  can  be  carried  out  for 
large  optical  depths.  Numerical  problems  that  may  occur  when  the  single¬ 
scattering  albedo  «  -  0  or  1  are  also  resolved  in  the  present  program.  Using 
results  from  the  adding  method  as  a  reference,  errors  in  the  5 -four- stream  scheme 
for  the  nonhomogeneous  multilayer  cases  remain  within  -5%.  Transfer  of  infrared 
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radiation  has  not  been  treated  satisfactorily  by  the  two-stream  techniques, 
because  they  do  not  proauce  reliable  results  in  the  limit  of  no  scattering.  This 
problem  is  circumvented  in  the  present  algorithm  by  using  the  "double  Gauss" 
quadrature  (Sykes,  1951),  which  produces  a  correct  relation  between  the  flux  and 
intensity  from  an  isotropic  source. 

Using  the  5-four-stream  scheme  along  with  the  parameterized  single¬ 
scattering  properties  reported  in  Section  2,  we  calculate  the  solar  reflectance 
and  transmittance  properties  for  typical  ice  clouds  with  different  thicknesses 
and  solar  zenith  angles.  These  properties  are  defined  by  r-Sr^w^  and  t-St^Wi 
where  r^  and  tj^  denote  the  reflectance  and  transmittance  (direct  plus  diffuse) 
for  an  individual  band,  and  w^  is  the  fraction  of  solar  flux  in  the  band.  The 
relative  errors  are  examined  with  respect  to  the  "exact"  results  computed  from 
the  adding  method  using  the  "exact"  phase  function.  Using  (Nov.  1)  as  an 
example,  Figs.  3a  and  3b  show  the  "exact"  reflectance  and  transmittance, 
respectively,  and  Figs.  3c  and  3d  show  the  corresponding  percentage  errors  in  the 
5 -four-stream  scheme  for  reflectance  and  transmittance.  The  5-four-stream 
results  in  most  cases  achieve  relative  accuracies  within  -5%.  The  errors  shown 
in  Figs.  3c  and  3d  contain  not  only  errors  due  to  the  5 -four- stream 
approximation,  but  also  errors  produced  by  the  parameterized  single-scattering 
properties  for  ice  clouds.  Generally,  errors  produced  by  the  latter  effect  are 
much  smaller  (<  1%). 

b.  Parameterizations  of  Nongray  Gaseous  Absorption 

In  each  band  shown  in  Table  2,  the  scattering  properties,  solar  irradiance, 
and  Planck  functions  may  be  treated  as  constant.  In  the  solar  spectrum, 
absorption  due  to  H2O  (2500-14500  cm"^)  ,  O3  (in  the  ultraviolet  and  visible),  CO2 
(2850-5250  cm”^),  and  O2  (A,  B,  and  7  bands)  is  accounted  for  in  the 
parameterization.  In  the  infrared  spectrum,  we  include  the  major  absorption 
bands  of  H2O  (0-2200  cm”^) ,  CO2  (540-800  cm'^) ,  O3  (980-1100  cm'^)  ,  CH*  (1100-1400 
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cm~^)  ,  and  N2O  (1100-1400  cm~^)  .  The  continuum  absorption  of  H2O  is  incorporated 
in  the  spectral  region  280-1250  cm”^. 

Nongray  gaseous  absorption  is  incorporated  into  the  6-four-stream  scheme 
based  on  the  correlated  k- distribution  method  developed  by  Fu  and  Liou  (1992) . 
In  this  method,  the  cumulative  probability,  g,  of  the  absorption  coefficient,  k^, 
in  a  spectral  interval  hv  is  used  to  replace  frequency,  u,  as  an  independent 
variable.  This  leads  to  an  immense  numerical  simplification,  in  which  about  ten 
thousand  frequency  intervals  can  be  replaced  by  a  few  g  intervals .  With  the 
required  accuracy  of  -1%,  we  find  that  the  minimum  niomber  of  quadrature  points 
in  the  g  space  ranges  from  one  to  about  ten  for  different  absorbing  gases  in 
different  spectral  regions.  In  order  to  treat  overlapping  of  CO2  and  H2O 
absorption,  the  correlated  k- distribution  has  been  formulated  for  a  single 
mixture  gas  to  economize  the  computational  time.  Using  a  minimum  number  of  g 
intervals  to  represent  the  gaseous  absorption  and  to  treat  overlap  within  each 
spectral  interval,  121  spectral  calculations  are  required  for  each  vertical 
profile. 

In  the  correlated  k-distribution  method,  the  vertical  nonhomogeneity  of  the 
atmosphere  is  accounted  for  by  assuming  a  simple  correlation  of  k-distribution 
at  different  temperatures  and  pressures.  The  accuracy  of  the  6-four-stream 
scheme  coupled  with  the  parameterization  program  for  nongray  gaseous  absorption 
has  been  examined  with  respect  to  results  computed  from  a  line-by-line  program. 
Shown  in  Fig.  4a  is  a  comparison  of  IR  cooling  rates  using  the  midlatitude  summer 
profile  when  all  the  constituents  are  included  in  the  calculations.  A  maximum 
heating  rate  difference  of  0.08  K/day  occurs  in  the  troposphere.  The  errors  in 
the  calculated  fl\ix  are  within  0.3%.  Figure  4b  shows  a  comparison  of  solar 
heating  rates  due  to  water  vapor  for  different  solar  zenith  angles.  The  errors 
in  heating  rates  are  less  than  0.05  K/day.  The  errors  jn  total  atmospheric 
absorption  and  downward  surface  fluxes  are  generally  less  than  0.5%. 
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c.  Molecular  Scattering 

The  average  extinction  coefficient  due  to  molecular  scattering  for  a  given 
band  may  be  expressed  as  (Slingo  and  Schrecker,  1982) 

^=Rp/T,  (3.2) 

where  p  is  pressure  (mb),  T  is  temperature  (K) ,  and  ^  is  in  units  of  m~^.  The 
coefficients  R  for  the  six  bands  in  the  solar  spectrum  are  9.022x10“®,  5.282x10"^, 
5.722x10“®,  1.433x10“®,  4.526x10“®,  and  1.529x10“®. 

For  Rayleigh  scattering,  the  single -scattering  albedo  is  1.  The  phase 
function  P(cos  6)  -  3/4(1  +  cos^  0),  so  that  the  expansion  coefficients  u>i  -  0, 

u>2  -  1/2,  and  w,  -  0,  where  i  >  3. 

d.  Single-Scattering  Properties  for  Combined  Cloud  Particles  and  Gases 
Consider  a  nonhomogeneous  atmosphere  which  is  divided  into  N  homogeneous 

layers.  For  the  case  when  both  Rayleigh  scattering  and  extinction  by  cloud 
particles  are  mixed  with  gaseous  absorption,  the  total  optical  depth  for  each 
layer  is 

4r(g)  =  At®  +  +  4r®(g)  ,  (3.3) 

where  Ar®  and  Ar”  represent  the  optical  depth  due  to  Rayleigh  molecules  and  cloud 
particles,  respectively.  Ar®(g)  is  the  optical  depth  contributed  by  the  gaseous 
absorption  for  a  given  g  (the  cumulative  probability) ,  which  can  be  expressed  by 

Ar®(g)  =  k(g,p,T)pAz,  (3.4) 

where  k(g,p,T)  is  the  equivalent  k-function  (Fu  and  Liou,  1992),  p  is  the  density 
of  the  absorber,  and  Az  is  the  geometric  thickness  of  the  layer.  In  this  case, 
Ar*’  -  At”  +  At”,  where  At”  and  At”  are  cloud  scattering  and  absorption  optical 
depths,  respectively.  Thus,  the  combined  single -scattering  albedo  and  the 
expansion  coefficients  for  the  phase  function  can  be  obtained  by 
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Ar“  +  Ar" 

Ar*^  +  Ar”  +  Ar^Cg) 


(3.5) 


Ar^w”  ■*•  AT*wf 


At,  +  Ar* 


(3.6) 


where  wj*  and  wf  denote  the  expansion  coefficients  for  the  phase  function  for 
cloud  particles  and  Rayleigh  molecules,  respectively.  Since  the  phase  function 
is  independent  of  gaseous  absorption,  the  combined  w,  is  constant  over  a  given 
spectral  absorption  band.  Once  the  single -scattering  properties  have  been 
defined  for  a  given  g  for  each  level,  monochromatic  radiative  transfer 
calculations  may  be  carried  out,  and  the  total  flux  over  each  spectral  band  can 
be  obtained  by  integrating  the  flux  solution  in  the  g  space.  As  described  in 
Section  3(b),  a  number  of  intervals  are  used  for  the  g-quadrature ,  which  are 
sufficient  to  confine  the  integration  errors  within  1%. 

4.  Comparisons  with  Observations 

Recently,  the  cirrus  intensive  field  observation  (IFO)  was  conducted  in 
Wisconsin,  from  13  October  to  2  November  1986,  as  a  major  component  of  the  First 
ISCCP  Regional  Experiment  (FIRE)  to  improve  our  understanding  of  cirrus  clouds 
and  their  influence  on  weather  and  climate  processes  (Starr,  1987).  In  cirrus 
IFO,  the  microphysical  measurements  were  simultaneously  obtained  along  with 
radiation  measurements  by  using  high  altitude  research  aircraft.  This  allows 
comparisons  to  be  made  between  model  calculations  and  in-situ  observations.  In 
this  section  we  present  comparisons  of  model  derived  infrared  emissivlty  and 
solar  albedo  with  those  obtained  from  the  NCAR  Sabreliner  on  28  October  1986. 


The  Sabreliner  was  flown  near  the  western  shore  of  Lake  Michigan  in  the 
vicinity  of  Green  Bay,  Wisconsin,  in  a  thin  banded  cirrus  layer,  between  1530  and 
1630  UTC  (0930-1030  local  time).  The  cloud  top  and  base  heights  were  about  11.1 
and  8.9  km,  respectively.  A  racetrack  pattern  was  flown  at  six  different  levels 
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throughout  the  cirrus.  The  mean  solar  zenith  angle  during  the  flight  was  61.3°. 
Since  radiative  and  microphysical  characteristics  of  the  cloud  sampled  on  the 
south  side  of  the  racetrack  (referred  to  as  cloud  1)  are  different  from  those 
sampled  on  the  north  side  (referred  to  as  cloud  2),  the  data  from  each  side  were 
analyzed  separately.  Furthermore,  each  side  of  the  racetrack  pattern  was  divided 
into  thinner,  mean,  and  thicker  cases  to  minimize  the  sampling  error.  Details 
of  the  description  of  the  observations  and  data  analysis  procedure  have  been 
given  by  Smith  et  al.  (1990),  and  Heymsfield  et  al.  (1990). 

The  effective  upward  emlssivity  and  downward  emissivity  £*  involving  a 
nonblack  cloud  may  be  defined  as  (Paltridge  and  Platt,  1981) 


and 


-  It  (Zb) 

C  ^  t 

oT*  -  Fj(Zb) 

£*  -  ^L(Zb)  -  Ft  (Zt) 
oT*  -  F^Uzt) 


(4.1) 


(4.2) 


where  F^  and  F^  are  the  infrared  broadband  upward  and  downward  fluxes,  and  z^ 
are  the  cloud  top  and  base  heights,  and  T  is  the  mean  cloud  temperature.  oT*  is 
the  flux  emitted  by  a  blackbody  at  T. 

The  solar  albedo  may  be  defined  as  (Paltridge  and  Platt,  1981) 

FsVzt)(l-a)  -  Fs*(Zb) 

Fs*(Zt)  -  Fs’(Zb)  (4.3) 


where 


^  ,  [Fs^(z,)-Fs\zt)]  -  [Fs*(Zb)  -  Fs^Zb)] 

Fs*(Zb)  +  Fs’(Zb) 


(4.4) 


In  these  equations,  Fs  is  the  solar  broadband  flux. 

e^,  «*,  and  a  were  derived  from  broadband  flux  observations  by  Stackhouse 
and  Stephens  (1991)  for  the  28  October  1986  FIRE  case,  which  are  shown  in  Figs. 
5a-c  as  a  function  of  the  ice  water  path  (IWP),  defined  as  IWC^Az.  The  error 


bars  through  these  observation  points  denote  the  uncertainties  of  the 
observation.  The  theoretical  results  are  displayed  in  these  figures  for 
comparisons.  In  the  calculations,  atmospheric  temperature  and  moisture  profiles 
were  averaged  over  the  two  soundings  from  Green  Bay  at  1400  and  1730  UTC.  For 
CO2,  CH4,  and  N2O,  uniform  mixings  were  assumed  throughout  the  atmosphere  with 
concentrations  of  330,  1.6,  and  0.28  ppmv,  respectively.  The  ozone  profile  and 
the  temperature  and  moisture  profiles  above  sounding  were  assumed  to  be  those  of 
the  US  Standard  Atmosphere.  The  cloud  was  positioned  between  8.9  and  11.1  km  and 
the  IWP  was  varied  by  changing  the  assumed  IWC  but  keeping  the  cloud  position 
fixed.  Ice  saturation  mixing  ratio  was  assumed  in  the  cloud.  In  the 
computations,  the  mean  effective  sizes  of  25,  50,  75,  100,  and  125  ^m  were  used, 
corresponding  to  five  solid  lines  for  c',  £*,  and  a  in  Figs.  5a-c.  For  solar 
radiation,  the  solar  zenith  angle  used  was  61.3*  and  surface  albedo  was  assumed 
to  be  0.072  (Stackhouse  and  Stephens,  1991).  Comparing  the  theoretical  results 
to  observed  data,  we  find  that  the  mean  effective  sizes  of  these  cases  lie  in  the 
range  of  -50  -  -100  nm  for  and  e*,  and  -50  -  -75  ^m  for  a.  The  mean  effective 
size  for  and  e*  in  each  case  appears  to  be  consistent  in  view  of  the 
observational  uncertainties.  However,  the  mean  effective  size  associated  with 
solar  albedo  is  smaller  than  that  associated  with  effective  emissivity.  This 
difference  can  be  explained  by  the  fact  that  solar  albedo  is  more  sensitive  to 
the  size  spectra  at  the  cloud  top.  Based  on  observations  (Heymsfield  et  al., 
1990) ,  ice  crystals  progressively  range  from  smaller  sizes  near  the  cloud  top  to 
larger  sizes  at  the  base  of  the  cloud.  For  the  same  cloud,  e.g.,  cloud  1,  Fig. 
5  shows  that  D,  Increases  as  IWP  Increases.  This  is  expected  because  larger  IWC 
is  usually  correlated  with  larger  particle  size.  As  given  by  Eqs.  (4.1)  and 
(4.2),  the  effective  emissivitles  are  defined  in  terms  of  fluxes  that  contain 
components  due  to  scattering  of  cloud  particles  and  absorption  of  gases, 
principally  water  vapor.  For  this  reason,  the  computed  effective  emissivities 
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may  exceed  1  when  IWP  is  large,  and  may  deviate  from  0  when  IWP  -*  0.  Figure  5d 
illustrates  the  plot  of  solar  albedo  vs.  effective  downward  emissivity  for  a 
number  of  mean  effective  sizes.  The  observed  data  derived  from  the  28  October 
1986  FIRE  case  are  also  displayed.  In  the  domain  of  albedo -emissivity,  the 
effect  of  ice  crystal  size  distribution  is  fairly  small.  Reasonable  agreements 
are  shown  between  the  theoretical  results  and  observations . 

We  also  calculate  the  emissivity  of  cloud  particles  alone  by  using  the  6- 
four- stream  scheme  with  the  parameterized  single -scattering  properties  described 
in  Section  2.  The  results  are  compared  with  those  derived  by  Smith  et  al.  (1990) 
for  the  28  October  1986  FIRE  case  in  Fig.  6,  where  each  pair  of  observed  points 
represents  the  emittance  from  the  top  layer  of  the  cloud  and  the  emittance  of  the 
entire  cloud  layer.  Smith  et  al.  utilized  a  broadband  infrared  radiative 
transfer  model  ^Gox  and  Griffith,  1979)  to  isolate  the  emissivity  of  ice 
particles  by  using  observed  downward  infrared  fluxes.  From  this  comparison  (Fig. 
6),  we  find  that  the  average  of  the  entire  cloud  layer  is  in  agreement  with 
the  value  derived  from  Fig.  5a  for  each  case.  Figure  6  also  shows  that  the 
averaged  mean  effective  size  in  the  upper  layer  of  the  cloud  is  generally  smaller 
than  that  of  the  entire  cloud  layer.  From  Fig.  6,  we  note  that  the  line  for  the 
D,  of  25  /im  crosses  the  50  line.  This  is  because  a  smaller  D,  produces  larger 
reflectivity,  r.  When  the  optical  depth  t  -♦  «,  the  emissivity  c-l-r  would  be 
smaller. 

From  Figs.  5  and  6,  the  estimated  D,  ranges  from  -50  to  -100  pim.  Using 
aspect  ratios  of  2  for  50  fim  and  3  for  100  /im,  the  effective  maximum  dimension 
(L,)  would  be  from  -100  to  -300  fjtm.  Wielicki  et  al.  (1990)  estimated  the  median 
mass  weighted  maximum  dimension  (L^,)  from  Sabreliner  observations  to  be  -75  to 
-400  ftm.  The  averaged  L^,  is  -238  ^m.  The  averaged  L,  from  the  present  study  is 
-200  ftm.  Note  that  and  L,  represent  volume  and  area  weighted  mean  values. 
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Nevertheless,  the  computed  L,  from  the  model  appears  to  compare  reasonably  well 
with  observed  data. 

In  light  of  preceding  discussion,  the  observed  cloud  emissivity  and  cloud 
albedo  can  be  reasonably  interpreted  from  the  theoretical  results ,  which  used  the 
scattering  and  absorption  properties  for  hexagonal  ice  crystals  (columns  and 
plates)  .  The  actual  ice  crystals  that  were  collected  by  the  optical  probe  aboard 
the  Sabreliner  during  the  FIRE  cirrus  experiment  were  largely  bullet  rosettes, 
columns,  plates,  and  compacted  spatial  particles  (occasionally  having  extensions) 
(Heymsfield  et  al.,  1990).  Although  the  exact  ice  crystal  shapes  were  not 
employed  in  the  scattering  and  absorption  calculations ,  the  theoretical  results 
account  for  the  general  hexagonal  structure  of  ice  particles  and  include  accurate 
fl\ix  calculations.  This,  perhaps,  explains  some  success  of  the  present 
parameterization  scheme  in  the  interpretation  of  the  observed  data  from  aircraft. 
Past  researchers  frequently  used  area- equivalent  spheres  to  interpret  the 
observed  radiative  properties  of  cirrus  (see,  e.g.,  Stackhouse  and  Stephens, 
1991) .  Since  equivalent  spheres  scatter  more  light  in  forward  directions  and 
have  smaller  single-scattering  albedo  than  nonspherical  ice  crystals,  the 
assumption  that  ice  particles  are  spheres  would  lead  to  an  underestimation  of  the 
solar  albedo  of  cirrus  clouds  (Kinne  and  Liou,  1989).  Increasing  the  amount  of 
small  spherical  ice  crystals  has  not  led  to  a  better  interpretation  of  the 
radiative  properties  of  cirrus  clouds.  Finally,  it  should  be  pointed  out  that 
further  verification  of  theoretical  results  would  require  narrowing  the  large 
experimental  uncertainties  that  are  related  to  radiation  and  microphysical 
measurements . 

Finally,  we  examine  the  relationship  for  the  emissivity  of  cirrus  clouds 
as  a  function  of  the  optical  depth  at  0.55  /im.  Using  the  present  model,  a  number 
of  emissivitles  were  computed  for  the  11  observed  cirrus  cloud  types  (Table  1) 
with  cloud  thicknesses  of  0.5,  1,  2  and  4  km.  The  results  are  plotted  in  Fig. 
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7  as  a  function  of  the  optical  depth.  The  solid  line  on  the  diagram  is  the  least 
square  curve  using  the  function 

€  =  1  -  exp(-  ar),  (4.5) 
where  r  is  the  visible  optical  depth.  The  fitted  coefficient  a  is  0.79.  Using 
a  diffusivity  factor  of  1.66,  the  visible  optical  depth  is  related  to  the 
infrared  optical  depth  (due  to  absorption)  by  a  factor  of  2.1.  This  result 
agrees  with  the  observed  result  reported  by  Platt  et  al.  (1987),  who  showed  that 
the  ratio  of  visible  extinction  to  infrared  absorption  varies  from  -1.8  to  2.3 
for  the  midlatitude  cirrus  clouds.  Note  that  the  11  size  distributions  used  in 
the  present  study  cover  a  spectrum  of  ice  crystal  size  spectra  that  occur  in 
midlatitude  cirrus  clouds. 

5.  The  Radiative  Effects  of  Ice  Crystal  Size  Distributions 

The  importance  of  cirrus  clouds  on  the  radiative  budget  of  the  earth  and 
the  atmosphere  has  been  documented  by  Liou  (1986).  Ackerman  et  al.  (1988)  have 
further  carried  out  an  investigation  regarding  the  internction  of  infrared  and 
solar  radiation  with  tropical  cirrus  anvils.  They  find  that  the  radiative 
effects  of  cirrus  clouds  are  very  sensitive  to  the  assumed  ice  water  content,  the 
cloud  thickness,  and  the  environmental  conditions  in  which  the  clouds  form. 
However,  their  study  was  confined  to  cirrus  clouds  with  a  fixed  size 
distribution.  In  this  study,  we  examine  the  effects  of  ice  crystal  size 
distributions  on  the  radiative  heating  and  budget  using  the  radiative  transfer 
model  developed  in  the  preceding  sections.  In  the  present  calculations,  standard 
atmospheric  temperature,  water  vapor  and  ozone  profiles  are  used.  The  cirrus 
cloud  layer  is  inserted  in  the  atmosphere  between  9  and  11  km.  For  water  vapor 
within  the  cloud  layer,  the  ice  saturation  mixing  ratio  is  assumed.  To  compute 
solar  fluxes,  a  solar  constant  of  1365  Wm~^,  a  surface  albedo  of  0.1,  a  solar 
zenith  angle  of  60",  and  a  12  h  solar  day  are  used.  The  computations  are  carried 
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out  by  using  the  mean  effective  sizes,  D,,  of  25,  50,  75,  100,  and  125  /im,  which 
cover  the  range  of  D,  from  observations, 
a.  Heating  Rate 

The  effects  of  the  ice  crystal  mean  effective  size  on  IR,  solar  and  net  (IR 
+  solar)  radiative  heating  rates  in  the  cloud  are  illustrated  in  Figs.  8a,  8b, 
and  8c,  respectively.  The  ItfP  used  is  30  g  m"^,  which  corresponds  to  an  IWC  of 
0.015  g  m"®.  In  the  calculations,  the  vertical  layer  thickness  is  set  to  be  250 
m  and  the  heating  rates  are  layer  averages.  For  the  smallest  D,  (25  fim) ,  IR 
heating/cooling  rates  of  about  30  K/day  take  place  close  to  the  bottom  and  top 
of  the  cloud  layer.  Such  large  differential  heating  rates  between  cloud  bottom 
and  top  (-60  K/day)  would  have  a  significant  effect  on  both  the  entrainment  and 
cloud  microphysics  (Ackerman  et  al.,  1988).  For  the  largest  D,  of  125  urn,  the 
IR  heating/cooling  rates  at  the  cloud  bottom  and  top  reduce  by  as  much  as  a 
factor  of  10.  A  cloud  consisting  of  smaller  D,  emits  more  IR  radiation  at  the 
cloud  top,  resulting  in  strong  cooling.  The  cloud  absorbs,  at  the  same  time, 
more  IR  radiation  emitted  from  the  warmer  surface  and  atmosphere  below,  resulting 
in  strong  heating  at  the  cloud  base  region.  If  IWC  is  fixed,  a  cloud  containing 
smaller  D,  would  have  larger  emittance/absorptance  because  of  larger  cross 
sectional  area.  The  solar  heating  rate  peaks  at  the  cloud  top  and  reaches  a 
value  of  -5  K/day  for  D,  of  25  nm.  The  net  heating  rates  due  to  both  IR  and 
solar  radiation  in  the  cloud  are  shown  in  Fig.  8c.  It  is  evident  that  the  cloud 
heating  rate  field  is  largely  controlled  by  IR  radiation  and  that  smaller  D,  can 
produce  larger  radiative  heating  gradients . 

For  a  given  ice  crystal  size  distribution,  larger  IWC  generates  a  stronger 
radiative  heating  gradient,  as  demonstrated  by  Ackerman  et  al.  (1988).  Based  on 
the  preceding  discussion  on  the  relation  of  D,  and  heating  rates,  changes  in  the 
heating  gradient  due  to  the  effects  of  D,  and  IWC  are  thus  in  opposite 
directions.  According  to  observations  (e.g.,  Heymsfleld  et  al .  ,  1990),  IWC 
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usually  (if  not  always)  Increases  with  increasing  ice  crystal  size.  Therefore, 
the  effects  of  IWC  and  D,  on  heating  rates  may  partially  be  compensated  for  if 
both  parameters  are  interactively  included  in  the  calculation, 
b.  Cloud  Forcing 

The  ice  crystal  size  distribution  of  cirrus  clouds  has  a  significant  effect 
on  the  radiation  budget  of  the  earth  and  the  atmosphere .  To  quantify  this 
effect,  we  use  the  cloud  radiative  forcing  concept  (e.g.,  Hartmann  et  al.,  1986) 
for  the  flux  at  the  top  of  the  atmosphere  (TOA)  and  define  the  following: 


^ir.« 


7  cl 
ir,» 


Fov 

ir,» ' 


where  ci  and  ov  denote  clear  and  overcast  conditions ,  and  and  F,  are  the 
upward  IR  and  solar  fluxes,  respectively.  The  net  cloud  radiative  forcing  is 


C  *  Cjr  +  C,, 

For  partly  cloudy  conditions,  if  the  cloud  cover  is  denoted  as  »>,  the  cloud 
radiative  forcing  is  given  by  i/C. 

Figures  9a- 9c  show  the  results  of  C,,  and  C  as  a  function  of  IWP  for 
different  mean  effective  sizes.  The  IR  cloud  forcing  (Fig.  9a)  is  always 
positive,  corresponding  to  the  heating  of  the  earth- atmosphere  system  due  to  the 
greenhouse  effect  of  clouds,  while  the  solar  cloud  forcing  C,  (Fig.  9b)  Is  always 
negative,  showing  the  cooling  of  the  system  by  the  cloud  albedo  effect. 

C,  are  strongly  dependent  on  both  IWP  and  D,.  For  an  IWP  of  30  g/m^,  C^i. 
decreases  from  131  to  5A  W/m^  as  D,  increases  from  25  to  125  ^m.  The  solar  cloud 
forcing  counterpart  increases  from  -114  to  -24  W/m^.  The  net  cloud  forcing  (Fig. 
9c)  shows  that  cirrus  clouds  have  a  net  heating  effect  on  the  system  except  for 
IWP  >  50  g  m~^  in  the  case  when  D,  -  25  ^m.  These  results  are  similar  to  those 
derived  from  a  climate  model  presented  by  Stephens  et  al.  (1990). 

Usually,  small  D,  is  associated  with  cold  cirrus  clouds  that  have  small 
IWP.  We  may  safely  assume  that  midlatitude  cirrus  clouds  would  have  positive 
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cloud  forcing,  i.e.,  the  greenhouse  effect  produced  by  the  presence  of  these 
clouds  is  more  pronounced  than  the  albedo  effect.  The  net  greenhouse  effect  is 
a  function  of  both  ice  crystal  size  and  IWP.  The  data  points  depicted  in  Fig. 
9c  are  computed  from  the  10  observed  ice  crystal  size  distributions  (except  Ci 
uncinus,  which  has  a  IWP  of  -220  g/m^)  using  a  thickness  of  2  km.  The  visible 
optical  depths  of  these  clouds  range  from  about  0.15  to  1.5.  Assuming  a  20% 
cloud  cover,  the  cloud  forcing  values  generally  range  from  about  4  to  8  W/m^. 
When  coupled  with  positive  radiative  forcings,  such  as  those  due  to  the  increase 
of  greenhouse  gases,  feedbacks  produced  by  the  increase  or  decrease  of  cirrus 
clouds  could  be  significant. 

Finally,  we  investigate  the  cloud  radiative  forcing  at  the  surface. 
Figures  10a  and  10b  show  the  IR  and  solar  cloud  forcing,  respectively.  In  the 
latter  case,  the  negative  cloud  forcing  values  produced  by  cloud  reflection  are 
substantially  similar  to  those  presented  in  Fig.  9b.  This  is  because  the 
atmosphere  is  largely  transparent  with  respect  to  solar  radiation  when  clouds  are 
present.  For  the  IR  case,  the  insertion  of  cirrus  clouds  produces  small  effects 
at  the  surface  (0-30  W/m^) ,  because  of  the  strong  absorption  of  the  atmosphere 
between  the  surface  and  the  cloud  in  the  infrared  wavelengths.  The  cirrus  cloud- 
induced  IR  heating  within  the  atmosphere  and  solar  cooling  at  the  surface  could 
have  a  significant  Impact  on  the  atmospheric  circulations  and  temperature 
distributions  (Ramanathan,  1987;  Sllngo  and  Slingo,  1988). 

6.  Summary  and  Conclusions 

In  this  paper,  a  new  approach  has  been  devised  for  the  parameterization  of 
radiative  transfer  in  the  atmosphere  involving  cirrus  clouds.  We  use  a  mean 
effective  size,  representing  ice  crystal  size  distribution,  in  the  development 
of  parameterization  for  the  single-scattering  properties  of  ice  crystals.  The 
extinction  coefficient,  the  single-scattering  albedo,  and  the  expansion 
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coefficients  of  the  phase  function  for  the  11  observed  ice  crystal  size 
distributions  have  been  computed  from  a  geometric  ray- tracing  technique  for 
hexagonal  ice  crystals  (size  parameter  >  30)  and  Mie-tjrpe  solution  for  spheroids 
(size  parameter  <  30).  Using  these  results  and  basic  physical  reasonings,  we 
develop  simple  polynomial  relationships  between  the  single -scattering  properties 
and  mean  effective  size  that  can  achieve  relative  accuracies  within  -1%. 

A  5 -four- stream  radiative  transfer  scheme  has  been  developed  for  flux 
calculations  in  both  solar  and  infrared  spectra.  With  respect  to  the  "exact" 
results  computed  from  the  adding  method  for  radiative  transfer,  this  scheme  can 
yield  relative  accuracies  within  -5%  for  all  atmospheric  conditions.  For  nongray 
gaseous  absorptions  due  to  H2O,  CO2,  O3,  CH4,  and  N2O,  we  use  the  correlated  k- 
distribution  method  to  obtain  a  number  of  equivalent  absorption  coefficients. 
A  total  of  121  spectral  calculations  for  each  vertical  profile  in  the  entire 
spectrum  is  required.  Compared  with  line-by-line  calculations  for  clear 
midlatitude  summer  atmosphere,  the  j-four-stream  scheme  coupled  with  the 
correlated  k-distribution  parameterization  produces  deviations  less  than  0.1 
K/day  in  heating  rates  and  within  0.5%  in  flux  calculations.  In  the  limit  of  no 
scattering,  the  5-four-stream  scheme  can  produce  reliable  results  since  the 
doubls  Gauss  quadrature  provides  a  correct  relation  between  the  flux  and 
intensity  from  an  Isotropic  source.  The  computational  time  of  the  scheme  is 
linear  with  respect  to  the  number  of  vertical  levels. 

Comparisons  have  been  carried  out  between  theoretical  calculations  and 
aircraft  measurements  involving  cloud  emissivlty  and  cloud  albedo  obtained  during 
the  cirrus  FIRE  IFO.  The  effect  of  ice  crystal  size  distribution  in  terms  of  the 
mean  effective  size  is  important  in  the  interpretation  of  the  observed  data. 
Further  verification  of  the  model  results  requires  the  reduction  of  large 
experimental  uncertainties  in  both  radiation  and  microphysical  measurements. 
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We  further  investigate  the  effects  of  ice  crystal  size  distributions  on 
the  radiative  properties  of  the  atmosphere  involving  cirrus  clouds.  For  a  given 
ice  water  path,  cirrus  clouds  containing  smaller  mean  effective  sizes  reflect 
more  solar  radiation,  trap  more  thermal  infrared,  and  produce  strong  cloud- top 
cooling  and  cloud-base  heating.  This  pattern  would  enhance  the  in- cloud  heating 
rate  gradient,  which  could  be  significant  in  terms  of  cloud  microphysical 
processes.  For  the  observed  midlatitude  cirrus  cases,  the  IR  greenhouse  effect 
outweighs  the  solar  albedo  effect.  The  degree  of  the  greenhouse  effect  involving 
cirrus  clouds  is  a  function  of  cloud  cover  as  well  as  the  mean  effective  size  and 
liquid  water  path.  Assuming  the  presence  within  a  sky  field  of  20%  cirrus 
clouds,  a  gain  of  -4-8  W/m^  is  obtained  from  the  present  calculation.  With 
respect  to  the  surface ,  the  presence  of  cirrus  clouds  reduces  the  solar  flux 
available  to  the  surface.  In  this  case,  the  solar  albedo  effect  dominates  the 
IR  greenhouse  effect,  which  is  confined  largely  within  the  atmosphere. 

There  are  a  number  of  significant  new  features  in  the  present  radiative 
scheme.  First,  the  ice  crystal  size  distribution  consisting  of  nonspherical 
particles  is  effectively  accounted  for  in  the  parameterization  of  the  single - 
scattering  properties  of  cirrus  clouds.  This  is  especially  important  in  the 
determination  of  the  solar  albedo  of  cirrus  clouds.  Second,  the  IR  and  solar 
radiation  regions  have  been  treated  in  a  self-consistent  fashion  using  the  same 
parameterization  scheme,  which  is  essential  to  the  investigation  of  the  solar 
albedo  versus  IR  greenhouse  effects  involving  clouds.  Third,  the  present  model 
can  efficiently  compute  th  ■.  detailed  vertical  structure  of  the  heating  rate 
profile  within  clouds,  which  is  critical  to  the  understanding  of  the  effects  of 
radiation  on  cloud  formation.  The  present  parameterization  for  the  radiative 
transfer  in  the  atmosphere  involving  cirrus  c?-ouds  is  well  suited  for 
incorporation  in  numerical  models  to  study  the  climatic  effects  of  cirrus  clouds. 


as  well  as  to  investigate  interactions  and  feedbacks  between  cloud  microphysics 
and  radiation. 
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Table  1.  Characteristics  of  the  11  ice  crystal  size  distributions  employed  in 


the  present  study. 

Particle  Size 
Distribution 

Number  of 
Particles  (cm“^) 

Ice  Water 
Content  (g  m"^) 

Mean  Effective 
Size  (/im) 

Cs 

0.187 

4.765  e-3 

41.5 

Ci  Uncinus 

0.213 

1.116  e-1 

123.6 

Ci  (Cold) 

0.176 

1.110  e-3 

23.9 

Ci  (Warm) 

0.442 

9.240  e-3 

47.6 

Ci  (T— 20'C) 

0.247 

8.613  e-3 

57.9 

Ci  (T— 40‘C) 

0.187 

9.177  e-3 

64.1 

Ci  (T— 60*0 

0.065 

6.598  e-4 

30.4 

Ci  (Oct.  22) 

0.052 

1.609  e-2 

104.1 

Ci  (Oct.  25) 

0.072 

2.923  e-2 

110.4 

Ci  (Nov.  1) 

0.032 

4.968  e-3 

75.1 

Ci  (Nov.  2) 

0.038 

1.406  e-2 

93.0 

Table  2.  Spectral  division  used  in  the  parameterization. 


Solar  Spectrum  Infrared  Spectrum 


Band  1 

Central  X 

(fim) 

Band  Limits 
(fim) 

Band  i 

Central  A 
(A»m) 

Band  Limits 
(cm“^) 

1 

0.55 

0.2-0. 7 

7 

4.9 

2200-1900 

2 

1.0 

0.7-1. 3 

8 

5.6 

1900-1700 

3 

1.6 

1.3-1. 9 

9 

6.5 

1700-1400 

4 

2.2 

1.9-2. 5 

10 

7.6 

1400-1250 

5 

3.0 

2. 5-3. 5 

11 

8.5 

1250-1100 

6 

3.7 

3. 5-4.0 

12 

9.6 

1100-980 

13 

11.3 

980-800 

14 

13.7 

800-670 

15 

16.6 

670-540 

16 

21.5 

540-400 

17 

30.0 

400-280 

280-10 


18 


70.0 


Table  3. 


Values  of  empirical  coefficients  in  Eq.  (2.9)*  for  the 
parameterization  of  the  extinction  coefficient  .  Note  that  the 
coefficients  are  determined  by  using  an  ice  density  of  0.9167  g 
cm"®.  For  other  ice  densities  p\,  they  should  be  adjusted  by  a 
factor  Pi/p\. 


Band  i 

ai 

^2 

1  -  6 

-6.656e-3 

3.686 

0.0 

7 

-7.770e-3 

3.734 

11.85 

8 

-8.088e-3 

3.717 

17.17 

9 

-8.441e-3 

3.715 

19.48 

10 

-9.061e-3 

3.741 

26.48 

11 

-9.609e-3 

3.768 

34.11 

12 

-1.153e-2 

4.109 

17.32 

13 

-8.294e-3 

3.925 

1.315 

14 

-1.026e-2 

4.105 

16.36 

15 

-1.151e-2 

4.182 

31.13 

16 

-1.704e-2 

4.830 

16.27 

17 

-1.741e-2 

5.541 

-58.42 

18 

-7.752e-3 

4.624 

-42.01 

*The  units  for  D, 

and  IVC  are  in  nm 

and  g  m"®. 

respectively. 

Table  4.  Values  of  empirical  coefficients  in  Eq.  (2.10)*  for  the 
parameterization  of  the  single-scattering  albedo  o). 


Band  1 

bo 

bi 

^2 

^3 

1 

.10998E-05 

-.26101E-07 

.10896E-08 

- .47387E-11 

2 

.20208E-04 

.96483E-05 

.83009E-07 

-  .32217E-09 

3 

.13590E-03 

.73453E-03 

.28281E-05 

- .18272E-07 

4 

-.16598E-02 

.20933E-02 

-.13977E-05 

-.18703E-07 

5 

.46180E+00 

.24471E-03 

-.27839E-05 

.10379E-07 

6 

.42362E-01 

.86425E-02 

-.75519E-04 

.24056E-06 

7 

.19960 

.37800E-02 

-.14910E-04 

.0 

8 

. 30140 

.26390E-02 

-.11160E-04 

.0 

9 

.39080 

.12720E-02 

-.55640E-05 

.0 

10 

.31050 

.26030E-02 

- .11390E-04 

.0 

11 

.20370 

.42470E-02 

-.18100E-04 

.0 

12 

.23070 

.38300E-02 

-.16160E-04 

.0 

13 

.56310 

- .14340E-02 

.62980E-05 

.0 

14 

.52070 

-.97780E-03 

.37250E-05 

.0 

15 

.32540 

.34340E-02 

-  .30810E-04 

.91430E-07 

16 

.10280 

.50190E-02 

-.20240E-04 

.0 

17 

.39640 

-.31550E-02 

.64170E-04 

- .29790E-06 

18 

.80790 

-.70040E-02 

.52090E-04 

-  .14250E-06 

♦The  units  for  D,  are  in  urn. 


Table  5a.  Values  of  empirical  coefficients  in  Eq.  (2.13)*  for  parameterization 
of  the  expansion  coefficients  for  the  phase  function  in  the  solar 
wavelengths . 


Band  1 

i 

^0,1 

^1,1 

^3,1 

1 

.22110E+01 

-.10398E-02 

.65199E-04 

- .34498E-06 

1 

2 

.32201E+01 

.94227E-03 

.80947E-04 

-.47428E-06 

3 

.41610E+01 

.74396E-03 

.82690E-04 

-.45251E-06 

4 

.51379E+01 

.51545E-02 

.11881E-04 

-.15556E-06 

1 

.22151E+01 

-.77982E-03 

.63750E-04 

-.34466E-06 

2 

2 

.31727E+01 

.15597E-02 

.82021E-04 

-.49665E-06 

3 

.40672E+01 

.25800E-02 

.71550E-04 

-.43051E-06 

4 

.49882E+01 

.86489E-02 

-.18318E-04 

- .59275E-07 

1 

.22376E+01 

.10293E-02 

.50842E-04 

-.30135E-06 

3 

2 

.31549E+01 

.47115E-02 

.70684E-04 

-.47622E-06 

3 

.39917E+01 

.82830E-02 

.53927E-04 

-.41778E-06 

4 

.48496E+01 

.15998E-01 

- .39320E-04 

-.43862E-07 

1 

.23012E+01 

.33854E-02 

.23528E-04 

. .20068E-06 

4 

2 

.31730E+01 

.93439E-02 

.36367E-04 

-.38390E-06 

3 

.39298E+01 

.16424E-01 

.10502E-04 

-.35086E-06 

4 

.47226E+01 

.25872E.01 

-.77542E-04 

-.21999E-07 

1 

.27975E+01 

.29741E-02 

-.32344E-04 

.11636E-06 

5 

2 

.43532E+01 

.11234E-01 

- .12081E-03 

.43435E-06 

3 

.56835E+01 

.24681E-01 

-.26480E-03 

.95314E-06 

4 

.68271E+01 

.42788E-01 

- .45615E-03 

.16368E-05 

1 

.19655E+01 

.20094E-01 

-.17067E-03 

.50806E-06 

6 

2 

.28803E+01 

.36091E-01 

-.28365E-03 

.79656E-06 

3 

.34613E+01 

.58525E-01 

- .46455E-03 

.13444E-05 

4 

.39568E+01 

.81480E-01 

-.64777E-03 

.19022E-05 

Band  i 

do 

di 

dj 

da 

1 

. 12495E+00 

-.43582E-03 

.14092E-04 

- .69565E-07 

2 

.12363E+00 

- .44419E-03 

.14038E-04 

- .68851E-07 

3 

.12117E+00 

-.48474E-03 

.12495E-04 

-.62411E-07 

4 

.11581E+00 

-.55031E-03 

.98776E-05 

-.50193E-07 

5 

- .15968E-03 

.10115E-04 

-.12472E-06 

.48667E-09 

6 

.13830E+00 

-.18921E-02 

.12030E-04 

- .31698E-07 

♦The  units  for  D,  are  in  pm. 


Table  5b.  Values  of  empirical  coefficients  in  Eq.  (2.15)*  for  parameterization 


of  the 

asymmetry  factor 

in  the  infrared 

wavelengths . 

Band 

i 

Co 

Cl 

C2 

c's 

7 

.79550 

2.524e-3 

-1.022e-5 

0.0 

8 

.86010 

1.599e-3 

-6.465e-6 

0.0 

9 

.89150 

1.060e-3 

-4.171e-6 

0.0 

10 

.87650 

1.198e-3 

-4.4856-6 

0.0 

11 

.88150 

9.858e-4 

-3.116e-6 

0.0 

12 

.91670 

5.499e-4 

-1.507e-6 

0.0 

13 

.90920 

9.295e-4 

-3.877e-6 

0.0 

14 

.84540 

1.429e-3 

-5.859e-6 

0.0 

15 

.76780 

2.571e-3 

-1.041e-5 

0.0 

16 

.72900 

2.132e-3 

-5.584e-6 

0.0 

17 

.70240 

4.561e-3 

-3.054e-5 

6.684e-8 

18 

.22920 

1.724e-2 

-1.573e-4 

4.995e-7 

♦The  units  for  D,  are  in  fim. 


Appendix 

The  5 -Four -Stream  Approximation  for  a  Homogeneous  Layer 


In  the  following,  we  present  the  necessary  equations  for  the  flux 
calculation  using  the  6 -four-stream  approximation.  The  notations  used  in  Liou 
et  al.  (1988)  are  followed.  Let  the  total  optical  depth  of  a  layer  be  r^.  The 
upward  and  downward  flttxes  at  a  given  optical  depth  r  (0<r<Ti)  are  given  by 

F*(r)  -  2ir(aiMili  +  a2fi2lz)  .  (A.l) 

F-(r)  - 


where  the  double  Gauss  quadratures  and  weights  in  the  four-stream  approximation 
are  -  -fi-i  —  0.2113248,  /12  “  “M-2  “  0.7886752,  and  a^  -  a_i  -  a2  “  a_2  “  0.5;  hq 
is  the  cosine  of  the  solar  zenith  angle;  and  xF©  denotes  the  direct  solar  flux 
at  the  top  of  the  atmosphere  (xFq  —  0  in  the  thermal  infrared) .  The  intensity 
solutions  can  be  written  in  the  form 


I2 

4>za2  ^lel  ^ie3  ^2®^ 

G2 

♦* 

Z2 

II 

I-l 

1-2 

- 

#2e2  #Jel  $ie3  $2®^ 
42®2  ^iel  ^^63  $2®^ 

^2e2  ^iel  ^JeS  ^2®^ 

Gi 

G-i 

G-2 

4- 

♦ 

Zl 

Zl 

Z2' 

(A.  3) 


where  fo(solar)  -  l/fiQ,  fo(thermal)  -  -l/riin(Bi/Bo) ,  Bq  and  Bj  are  Planck 
functions  evaluated  at  the  top  and  bottom  of  the  layer,  respectively,  and  the 
other  terms  except  2  defined  as  follows: 


el  -  e2  -  e'‘2\  e3  *  e4  * 

Zl,2  ”  j(Vl,Z  -  *ll,z)  1 

.(1± 

^  a 

<2  =  4(Ai,2  ±  2) . 

where  ki_2.  A^  2.  Vi.Zf  i,z>  ^22 »  ^12  calculated  from  the 


following  equations,  which  are  in  the  order  of  calculation.  The  single- 


scattering  albedo  is  denoted  by  w,  and  W|  (i”0, 1 , 2 , 3)  is  the  expansion  coefficient 
of  the  scattering  phase  function: 

3 

bi(solar)  -  “  jrFoJ^  (Mi)P,  (-po)//*! 

« -0 

b^( thermal)  =  (1  -  i  =  -2, -1,1, 2, 

Cj^j  “  (Mj)  »  1-iJ  “  “2  , -1 , 1 , 2 


b,‘ 


f  Cij/Pi, 

i*- j 

|(C2j  -  1)/P2, 

i  =  j. 

• 

’  b22  -  b2-2  b|2 

*>  b22  ~  b2.2 

'  b22  -  b2-2  bf2 

*  b22  i  bj.2 

• 

b2  ±  b.2  bf 

-  b2  ±  b.2. 

*  b22b22  b22b2i 

®21  “  b22b2l  +  b22b2i 

■  b22b22  b2ib22 

®  b22b22  +  b2ib2i 

• 

d2  •  b22b2  +  b2ibi  +  b2fo 
“  ^12^2  +  bi^bi  ■*■  bj^fg, 

®22  *  ^22^22  ^12^21  ^21  ~  ^22^21  *  b2ibii 

®12  *  bi2b22  *  biib22  an  =  bi2b22  +  b2ib2i 

^2  ”  b22b2  +  b22b2  +  bj  f  g 
di  *  b22b2  +  b22b2  +  b2fg. 


b  ■  a22  a22  c  *  322322  ~  322322  3  =  b22b2i  ~  b22b2i, 

k2  *  [(b+\A>2+4c  )/2]^'2 

^1,2  *  (^1,2  “  322)/322» 


k2  =  [(b-v^2+4c  )/2]1/2 


*/l  ®  (*^1^0  ■*■  ®12^2  ~  *22^1  )/^ 


»?2  *  (^2^0  +  ^21^1  -  aiid2)/f' 
12  ~  (^2^0  ■'■  ^21*^1  ~  ^11^2)/^ 


f'  -  fo*  -  bfo^  -  c. 


The  coefficients  0^(1  -  ±1,2)  are  to  be  determined  from  radiation  boundary 
conditions.  Consider  a  homogeneous  cloud  layer  characterized  by  an  optical  depth 
Ti,  and  assume  that  there  is  no  diffuse  radiation  from  the  top  and  bottom  of  this 
layer,  then  the  boundary  conditions  are 

I-i,-2<*’  =  0)  -  0^ 

h.2('  =  ^1)  =  0  r  (A. 4) 


The  boundary  conditions  can  be  modified  to  include  the  nonzero  diffuse  radiation. 
Gi  can  be  obtained  by  an  inversion  of  a  four-by-four  matrix  in  Eq,  (A. 3).  We 
note  that  the  only  difference  in  the  four- stream  formulation  between  infrared  and 
solar  wavelengths  is  the  definition  of  fo  and  b^  (1  '•  ±1,2). 

It  is  possible  to  incorporate  a  5-function  adjustment  to  account  for  the 
forward  diffraction  peak  in  the  context  of  the  four-stream  approximation.  We  may 
use  the  similarity  principle  for  radiative  transfer  to  adjust  the  optical  depth, 
single-scattering  albedo  and  expansion  coefficients  of  the  phase  function  in  the 
forms  (Liou  et  al.,  1988) 

r'  •  r(l  -  fw) 

«  *  (1  -  f)w/(l  -  fw) 

wi  -  [w,  -  f(2i  +!)]/(!  -  f),  1,2,3, 


where  the  fraction  of  scattered  energy  residing  in  the  forward  peak,  f  -  u^/9 . 


Figure  Captions 


Fig.  1  (a)  Extinction  coefficient/IWC,  (b)  single -scattering  albedo  and  (c) 

asymmetry  factor  as  functions  of  the  mean  effective  size  for  the 
spectral  interval  1.9-2. 5  pm.  The  "0"  points  represent  exact 
results  from  the  light-scattering  program,  and  the  curves  are  from 
the  parameterizations . 

Fig.  2  Same  as  Fig.  1,  except  for  the  spectral  interval  800-980  cm"^. 

Fig.  3  (a)  Solar  reflectance  and  (b)  solar  transmittance  calculated  by  the 


adding  method 

using  "exact”  phase 

function  for 

Ci 

(Nov.  1). 

The 

corresponding 

percentage 

errors 
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reflectance 

and 

transmittance 
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scheme 

with 

parameterizations  of  the  single -scattering  properties  are  displayed 
in  (c)  and  (d) . 

Fig.  4  Comparisons  between  the  present  model  and  line-by-line  calculations 

for  (a)  infrared  cooling  rates  due  to  H2O  (with  continuum)  +  CO2  + 
O3  +  CH*  +  N2O  in  the  spectral  region  10-2200  cm~^,  and  (b)  solar 
heating  rates  due  to  H2O  in  the  spectral  region  2500-14500  cm"^,  in 
the  midlatitude  summer  atmosphere. 

Fig.  5  (a)  Effective  upward  emissivity,  (b)  effective  downward  emissivity, 

and  (c)  solar  albedo  as  functions  of  the  ice  water  path  (IWP). 
Model  results  representing  the  mean  effective  sizes  of  25,  50,  75, 
100  and  125  pm  are  compared  with  the  observed  values  obtained  from 
the  cirrus  FIRE  IFO. 

Fig.  5(d)  Solar  albedo  as  a  function  of  effective  downward  emissivity.  Model 
results  for  the  mean  effective  sizes  of  25,  50,  75,  100,  and  125  pm 
are  compared  with  the  observed  values  obtained  from  the  cirrus  FIRE 
IFO. 
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Fig.  6 


Fig.  7 


Fig.  8 


Fig.  9 


Fig.  10 


The  emissivity  of  cloud  particles  alone  as  a  function  of  the  ice 
water  path  (IWP).  Model  results  for  the  mean  effective  sizes  of  25, 
50,  75,  100  and  125  /tm  are  compared  with  the  observed  values 

obtained  from  the  cirrus  FIRE  IFO. 

The  emissivity  of  cirrus  clouds  as  a  function  of  the  visible  optical 
depth.  The  "x"  points  are  obtained  from  the  present  model  for  the 
11  cloud  types  (Table  1)  with  thicknesses  of  0.5,  1,  2  and  4  km. 
The  solid  line  is  from  least-square  fits. 

Heating  rates  in  a  cirrus  cloud  with  an  ice  water  content  of  0.015 

g/m^  for  five  mean  effective  sizes  of  25  (-  -),  50  ( - ),  75  ( - ), 

100  ( - ),  and  125  fim  ( - ).  (a)  Infrared,  (b)  solar,  and  (c) 

net.  For  solar  radiation,  the  solar  zenith  angle  is  60°,  the 
surface  albedo  is  0.1,  and  a  12  h  solar  day  is  used. 

Cloud  radiative  forcing  at  the  top  of  the  atmosphere  as  a  function 
of  the  ice  water  path  for  five  mean  effective  sizes  of  25  (-  -),  50 

(--),  75  ( - ),  100  ( - ),  and  125  /sm  ( - ).  (a)  Infrared,  (b) 

solar,  and  (c)  net.  For  solar  radiation,  the  solar  zenith  angle  is 
60°,  the  surface  albedo  is  0.1,  and  a  12  h  solar  day  is  used. 

Same  as  Fig.  9,  except  for  the  cloud  radiative  forcing  at  the 
surface,  (a)  Infrared  and  (b)  solar. 
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